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CHAPTER  1 
INTRODUCTION 

1,1  Motivation 

Design  of  mechanisms  and  machines  is  of  major  importance  in  mechan¬ 
ical  design.  Design  techniques  currently  in  use  are  generally  oriented 
toward  design  of  specific  small-scale  mechanisms.  Some  assumptions  made 
in  these  techniques  are  often  unrealistic.  The  increasing  complexity  of 
mechanisms,  as  required  by  machines  with  automatic  and  programmable 
action,  requires  general-purpose  techniques  for  design  of  large-scale, 
multidegree-of-f reedom  mechanisms  * 

The  purpose  of  this  research  is  to  develop  a  general-purpose  theory 
for  design  optimization  of  mechanisms  and  machines.  An  experimental 
computer  code  to  implement  the  theory  is  developed  and  a  number  of 
example  problems  are  considered  to  demonstrate  the  wide  range  of 
applicability  of  the  technique. 

1,2  Existing  Methods  of  Mechanism  Synthesis 
Synthesis  techniques  for  mechanisms  have  generally  been  aimed  at 
designing  systems,  some  member  of  which  is  required  to  either  describe  a 
desired  path  or  generate  a  function  of  the  input.  The  objective  in 
these  design  situations  is  to  determine  the  member  lengths  and  other 
mechanism  parameters  to  minimize  the  difference  between  the  desired  and 
actual  path  or  function  generated  by  the  mechanism.  The  precision  point 


approach  and  its  variants  [1,2]  have  been  applied,  with  success,  for 
solution  of  these  design  problems.  The  basic  idea  underlying  these 
approaches  is  to  design  the  mechanism  to  ensure  conformity,  within  cer¬ 
tain  tolerances,  between  the  actual  path  and  the  desired  path,  at  speci¬ 
fied  points  on  the  path.  The  Chebyshev  and  least  square  error  criteria 
have  been  used  to  characterize  the  error. 

Balancing  of  mechanisms  and  machines  is  another  area  that  has  re¬ 
ceived  considerable  attention  [3,4,5].  The  mechanisms  being  considered 
in  these  investigations  are  generally  high-speed  inertia  variant  rota¬ 
ting  machinery.  Due  to  the  inertia  variant  nature  of  the  mechanism,  the 
support  frame  experiences  large  shaking  forces  and  moments.  The  design 
objective  is  to  redistribute  mass  of  links  or  to  add  counterweights  to 
minimize  shaking  forces  or  moments. 

Synthesis  of  multidegree-of-f reedom  mechanisms  has  been  the  sub¬ 
ject  of  several  papers  [6,7].  The  synthesis  technique  described  in 
these  references  was  essentially  restricted  to  two-degree-of -freedom 
mechanisms.  Synthesis  of  a  mechanism  so  that  it  will  occupy  less  than 
a  prescribed  amount  of  space  has  been  studied  in  Ref.  8. 

In  the  synthesis  methods  considered  above,  constraints  have  gener¬ 
ally  been  imposed  on  design  variables,  such  as  link  lengths,  or  on  gen¬ 
eralized  coordinates,  such  as  angles  between  links.  Constraints  on 
transmission  angle  have  been  extensively  used.  Methods  for  stress  and 
deformation  constrained  design  of  mechanisms  have  recently  appeared  in 
the  literature  [9].  Minimum  weight  design  is  the  objective  of  these 
design  schemes. 


As  is  evident  from  the  brief  survey  given  in  the  preceding  para¬ 
graphs,  most  available  synthesis  schemes  are  oriented  toward  design 
of  a  specific  type  of  mechanism,  to  perform  a  specific  class  of  tasks. 
Some  efforts  have  been  made  in  developing  synthesis  schemes  that  are 
more  general  than  those  described  above  [10,11,12].  The  generality  of 
these  schemes,  however,  is  limited  to  a  particular  class  of  problems. 

1.3  Modeling  Techniques  for  Large-Scale 
Mechanisms  and  Machines 

Modeling  techniques  for  large-scale  dynamic  mechanical  systems  have 
been  developed  only  in  the  1970’s  [13].  Modeling  techniques  for  dynamic 
electronic  and  structural  systems,  on  the  other  hand,  have  been  avail¬ 
able  for  some  time.  Two  of  the  modeling  methods  for  dynamic  mechanical 
systems  are  considered  appropriate  for  modeling  kinematic  mechanical 
systems.  One  of  them,  the  loop  closure  method  [14],  is  embodied  in  the 
computer  code  IMP  [15].  This  modeling  method  has  been  commonly  used  for 
modeling  kinematic  systems  [16]*  The  other  method,  the  ’constrained 
multielement  formulation’,  is  the  basis  for  computer  codes  ADAMS  [17] 
and  DADS  [18,19].  This  modeling  method  involves  writing  equations  of 
motion  for  each  individual  member  and  then  adjoining  equations  of  con¬ 
straint  through  Lagrange  multipliers.  This  modeling  method,  though  not 
yet  used  for  synthesis  of  kinematic  systems,  has  attractive  features  for 
doing  so. 


1.4  Techniques  Available  for  Design  Optimization 
of  Large-Scale  Systems 

Techniques  for  design  optimization  play  an  important  role  in 
kinematic  design.  Most  methods  previously  employed  for  optimization  of 
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structural  and  mechanical  systems  belong  to  the  field  of  nonlinear 
programming,  in  which  the  design  problem  is  formulated  in  terms  of 
design  variables  that  are  to  be  selected*  Optimization  methods  such  as 
the  Sequential  Unconstrained  Minimization  Technique  (SUMT)  [8]  and 
Optimality  criteria  [9]  have  been  used  for  Kinematic  synthesis. 
Performance  constraints,  however,  are  most  naturally  stated  in  terms  of 
state  or  response  variables.  Ad  hoc  techniques  have  been  used  to  reduce 
the  design  problem  to  a  standard  nonlinear  programming  problem,  with 
attendant  limitations  on  generality,  analytical  feasibility,  and 
computational  efficiency. 

Numerical  methods  used  in  optimal  control  and  optimal  design  theory 
[Ref.  20,  Ch.  3]  sharply  contrast  those  employed  in  early  mechanical 
system  optimization.  A  state  space  formulation  is  employed  that  explic¬ 
itly  treats  design  and  state  variables.  The  state  variable  is  generally 
the  solution  of  a  matrix  or  differential  equation,  for  which  an  adjoint 
or  costate  variable  is  defined  as  the  solution  of  a  related  problem. 

The  adjoint  variable  and  the  associated  adjoint  equation  are  used  to 
provide  explicit  design  sensitivity  information.  This  sensitivity 
information  is  required  for  virtually  all  iterative  methods  of  design 
optimization.  The  state  space  optimization  technique  has  been  success¬ 
fully  applied  to  design  of  structural  systems.  When  applied  to  optimi¬ 
zation  of  large-scale  structural  systems,  this  method  compares  favorably 
to  indirect  methods  such  as  optimality  criteria  [20]. 

Most  mechanical  design  problems  require  the  system  to  perform  over 
the  range  of  input  or  control  parameters.  In  the  case  of  kinematic  syn¬ 
thesis  of  a  four-bar  path  generator  mechanism,  for  example,  one  could 
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consider  the  angular  input  given  to  the  input  link  to  be  the  input 
parameter.  At  every  point  in  the  specified  range  of  the  input  parameter, 
design  constraints  must  be  satisfied.  Most  kinematic  synthesis  problems 
fall  into  this  class,  called  worst  case  or  parametric  optimal  design. 

This  optimal  design  scheme  has  been  applied  with  success  to  design  of 
structural  systems  and  vehicle  suspension  systems.  [Ref.  20,  Ch.  5]. 

1.5  Scope  of  The  Report 

In  light  of  the  comments  made  in  Sections  1.1  to  1.4,  the  general 
conclusion  can  be  drawn  that  no  general-purpose  techniques  have  been 
developed  for  mechanism  synthesis  and  design.  The  level  of  generality 
implied  in  the  term  ’general  purpose1  is  the  ability  of  techniques  to 
handle  large-scale  mechanisms,  with  a  variety  of  constraints  and  cost 
functions  imposed  on  the  design. 

In  this  report,  a  technique  based  on  the  constrained  multielement 
formulation  is  developed  for  modeling  planar  kinematic  systems.  The 
technique  is  general  and  is  capable  of  modeling  multidegree-of-freedom 
systems.  Velocities  and  accelerations  of  the  members  and  reaction 
forces  in  joints  can  also  be  computed  and  constrained.  A  kinetostatic 
force  analysis,  as  opposed  to  a  time  response  analysis,  is  resorted  to 
for  computation  of  the  reaction  forces.  Since  the  basic  assumption  of 
the  kinetostatic  analysis  is  that  the  kinematics  of  the  system  are 
independent  of  externally  applied  forces,  the  dynamic  effects  due  to 
inertia  of  the  members  are  also  included  in  the  kinetostatic  force 


analysis 


The  state  space  optimization  technique  is  used  to  develop  a  general 
method  for  design  sensitivity  analysis.  The  method  allows  constraints 
to  be  imposed  on  functions  of  design  and  state  variables.  The  design 
sensitivity  information  so  obtained  is  used  in  the  gradient  projection 
method  for  iterative  optimization. 
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CHAPTER  2 

KINEMATIC  ANALYSIS  OF  MECHANISMS 

2.1  Introduction  to  the  Constrained 
Multielement  Formulation  ~~ 

Before  any  mechanism  synthesis  and  design  schemes  are  considered,  a 
technique  for  kinematic  analysis  of  mechanisms  needs  to  be  developed. 
Development  of  a  general  synthesis  and  design  procedure  requires  that 
the  analysis  procedure  be  general  in  nature*  Techniques  used  to  date 
have  been  based  on  writing  algebraic  equations  for  independent  loops  in 
a  mechanism.  These  equations  typically  involve  relative  position  var¬ 
iables  of  the  links  of  the  mechanism.  These  techniques ,  though  adequate 
for  analysis  of  closed— loop  mechanisms ,  are  not  well  suited  for  analysis 
of  open-loop  mechanisms. 

A  technique  that  has  been  very  successful  in  modeling  large-scale 
open—  and  closed— loop  mechanical  systems  is  the  constrained  multielement 
method  [17,18].  This  technique,  though  not  yet  used  for  kinematic 
modeling  of  mechanisms,  has  attractive  features  for  doing  so.  A  method 
for  kinematic  analysis  of  mechanisms  based  on  this  technique  is  devel¬ 
oped  in  the  following  sections. 

The  modeling  philosophy  of  the  constrained  multielement  (CME)  for¬ 
mulation  is  to  embed  a  local  coordinate  system  in  each  link  of  the  mech¬ 
anism  or  machine.  The  location  of  this  coordinate  system  is  arbitrary, 
but  is  often  located  at  the  center  of  mass.  Since  only  planar  systems 
are  being  considered  in  this  report,  the  position  and  orientation 
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of  any  body  or  link  in  the  system  can  be  described  by  the  three  gener¬ 
alized  coordinates  y^  and  Oi.  These  quantities  can  be  represented 
by  the  subvector  q(i),  where 

IYI 


(i) 

q 


(2-1) 


As  shown  in  Fig.  2.1,  any  point  p  on  body  i  of  the  system  can  be 
represented  by  coordinates  £  and  ri  ,  conveniently  measured  in  the 

ij  ij 

local  coordinate  system.  Knowing  the  generalized  coordinates  for  the 
body ,  it  is  then  possible  to  express  the  position  of  a  typical  point  p 
in  the  global  coordinate  system. 

A  mechanical  system  generally  consists  of  many  members,  connected 
by  joints.  These  joints  could  be  looked  upon  as  constraints  imposed  on 
the  relative  motion  of  connected  pairs  of  bodies.  The  CME  formulation 
thus  represents  joints  as  constraints  between  bodies  making  up  the 
system.  These  constraints  are  expressed  as  algebraic  equations  invol¬ 
ving  the  generalized  coordinates  of  the  two  connected  bodies  and  addi¬ 
tional  geometric  variables. 


2.2  Position  Analysis  Of  Mechanisms 
2.2.1  Formulation  of  State  Equation  for  Position 
2.2.1. 1  Kinematic  Equations  of  Constraint 


* 


* 


Consider  a  general  system  of  n  bodies  connected  by  i  independent 
joints.  Two  types  of  joints  (revolute  and  translation)  are  considered 
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in  this  report,  each  giving  rise  to  two  equations  of  constraint.  These 
equations  of  constraint  express  conditions  that  the  motion  of  the  two 
connected  bodies  must  satisfy  to  be  compatible  with  the  joint. 

A  system  of  n  bodies  in  the  plane  has  a  total  of  3n  generalized 
coordinates.  However,  if  this  system  has  £  independent  joints,  there 
are  2£  equations  of  constraint  between  the  3n  generalized  coordinates. 
Thus  the  number  of  f ree-degrees-of-f reedom  can  be  written  as 

m  =  3n  -  2£  (2.2) 

where 


n  =  number  of  bodies  in  the  system 

£  =  number  of  independent  joints  in  the  system 

The  condition  m  >  0  must  be  satisfied  by  all  kinematic  systems. 

The  generalized  coordinates  for  the  entire  system  can  be  denoted  by  the 

3n 

position  state  vector  z  e  R  ,  where 


z 


(2.3) 


Assuming  that  the  system  has  a  set  of  specified  design  variables 
s 

b  e  R  ,  the  kinematic  equations  of  constraint  can  be  written  as 


k 

$  (z,b)  =  0 


where 


(2.4) 
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k 

$  (z,b)  = 


<f^(z,b) 


k 

<(»  (z,b) 
2 


k 

<f>  (z,b) 

21 


<b  (z,b)  is  the  first  kinematic  constraint  due  to  joint  i  and  <\>  (z,b) 

2i-l  2i 

is  the  second  kinematic  constraint  due  to  joint  i. 


2 .2 .1.2  Kinematic  Driving  Equations 

Equation  2.4  as  presented  in  the  preceding  section,  is  a  system  of 
2%  equations  in  3n  variables.  Since  for  a  kinematic  system  3n  >  2ly  Eq. 
2.4  is  a  system  of  fewer  equations  than  unknowns.  To  solve  for  z  from 
Eq.  2.4,  m  additional  equations  are  required.  These  equations  can  be 
developed  by  observing  that  the  basic  purpose  of  designing  kinematic 
systems  is  to  obtain  a  system  that  transmits  motion  from  input  links  to 
output  links.  The  mechanism  or  machine  can  only  be  given  input  motion 
through  a  set  of  free  degrees  of  freedom.  These  free  degrees  of  freedom 
can  thus  be  specified  as  functions  of  some  free  parameter,  or  may  be 
specified  by  some  relationship  between  the  3n  state  variables.  Since 
all  the  free  degrees  of  freedom  must  be  specified,  to  drive  the  mechan¬ 
ism  in  a  unique  way,  m  additional  driving  equations  of  constraints 
arise,  in  the  form 


d 

$  (z,b,ot)  =  0 


(2.5) 
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where 


d 

$  (z  ,b ,  a) 


d 

<|^(z,b,a) 


d 

<{>  (z  ,b ,  a) 
m 

P  d 

a  e  R  is  a  vector  of  input  parameters,  and  cj)  (z,b,a)  represents  the 

i 

ith  driving  equation. 


2. 2. 1.3  State  Equation  for  Position 

Combining  Eqs.  2.4  and  2.5,  one  has  3n  independent  equations,  which 
may  be  written  as 


$>  (  z  ,b ,  a)  = 


k 

$  (z,b) 


=  0 


(2.6) 


_J>  (z,b,a)J 

Equation  2.6  is  the  state  equation  for  position  of  the  mechanism. 

Specifying  the  design  variable  vector  b  and  the  input  parameter  vector  a 

makes  Eq .  2.6  a  system  of  3n  independent  equations  in  3n  unknowns, 

3n 

z  e  R  .  Since  these  equations  are  highly  nonlinear,  more  than  one 
solution  for  z  is  possible.  Conversely,  for  some  designs  and  inputs, 
no  solution  may  exist. 


2.2.2.  Kinematic  Constraint  Equations 
Before  any  technique  is  considered  for  the  solution  of  Eq.  2.6,  it 
is  necessary  to  determine  the  explicit  form  of  the  equations  constitu¬ 
ting  this  system  of  equations.  Since  kinematic  equations  of  constraint 
occur  in  a  general  form;  i.e.,  the  equations  of  constraint  for  all 


joints  of  the  same  type  have  the  same  general  form,  it  is  sufficient  to 
consider  a  typical  joint  of  each  type.  Since  this  report  considers  only 
revolute  and  translation  joints,  a  typical  joint  of  each  type  will  be 
treated.  Driving  constraints,  on  the  other  hand,  do  not  lend  them¬ 
selves  to  the  same  general  characterization.  These  constraints  tend 
to  be  problem  dependent. 


2.2.2. 1  Constraint  Equations  for  a  Revolute  Joint 

Figure  2.2  shows  the  two  adjacent  bodies  i  and  j,  with  body-fixed 
coordinate  systems  0  x  y  and  0  x  y  ,  respectively.  The  origins  of 

i  i  i  j  j  j 

these  reference  frames  are  located  in  the  global  reference  frame  by 
vectors  R  and  R  ,  respectively.  Let  point  p  on  body  i  be  located  by 

1  j 

a  body-fixed  vector  r  and  point  p  on  body  j  be  located  by  a  body 

ij  ji 

fixed  vector  7  .  Points  p  and  p  are,  in  turn,  connected  by  a 

ji  ij  ji 

vector  r  .  Executing  a  closed  path  from  the  origin  of  the  fixed 
P 

reference  frame  yields  the  vector  relationship. 

R  +  7  +7  -7  -  R  =  O’  (2.7) 

i  ij  P  ji  j 

Demanding  that  points  p  and  p  are  coincident  guarantees  the 

ij  ji 

existence  of  a  rotational  joint  between  bodies  i  and  j  at  this  common 

point.  This  is  equivalent  to  the  condition  r  =0.  Thus,  Eq.  2.7 

P 

becomes 


R+r  -r  -  R  =  0  (2.8) 

i  ij  ji  j 


In  matrix  form,  Eq.  2.8  can  be  written  as 
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is  a  rotation  matrix  from  the  local  to  global  reference  frame. 

Expanding  Eq .  2.9,  the  equations  of  constraint  for  the  revolute 
joint  can  be  written  as 

A  =  x  +  £  cos  0  -  n  sin  0  -  x  -  £  cos  0  +  n  „  sin  0=0 

x  i  ij  i  ij  i  j  ji  j  J1  3 

<f>  =  y  +  £  sin  0+q  cos  0  -  y  -  ?  sin  0  -  n  cos  0  -  0 

y  i  ij  i  ij  i  j  ji  j  ji  J 

(2.10) 

In  Eq.  2.10,  x  ,  y  ,  0  ,  x  ,  y  ,  and  0  are  state  variables.  The 

i  i  i  j  j  j 

variables  £  ,  n  ,  £  ,  and  n  are  related  to  the  length  of  the 

ij  ij  ji  ji 

members  and  hence  to  the  design  variables. 

2. 2. 2. 2  Constraint  Equations  for  a  Translational  Joint 

Figure  2.3  shows  two  bodies  connected  by  a  translational  joint. 

For  this  type  of  joint,  the  points  p  and  p  lie  on  a  line  parallel  to 

ij  ji 

the  path  of  relative  motion  between  the  two  bodies.  These  points  are 

located  by  nonzero  body-fixed  vectors,  r  and  r  ,  that  are  perpen- 

ij  Ji 

dicular  to  the  line  of  relative  motion.  A  scalar  equation  of  constraint 

can  be  written  by  taking  the  dot  product  of  r  with  r  .  Since  these 

ij  P 

two  vectors  are  perpendicular,  their  dot  product  must  vanish. 
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^Figure  2.3  Translation  Joint 
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r  •  r  =0  (2.11) 

P  ij 


Using  Eq.  2.7  to  solve  for  r  ,  Eq.  2.11  becomes 

P 


(R+r  )  -  (R  +  r  )  •  r  =  0  (2.12) 

J  ji  i  ij  ij 


Using  the  transformation  matrix  to  express  r  and  r  in  the  global 

ji  ij 

reference  frame,  Eq .  2.12  can  be  written  as 


=  (u 

-  X 

)(U 

- 

U 

)  +  (V 

> 

V-/ 

/*“\ 

1 

-  V  )  =  0 

(2.13) 

n 

i 

i  i 

j 

i 

i  i 

j 

where 

U 

-  X 

+ 

5 

cos 

6 

_ 

T1 

sin 

e 

i 

i 

ij 

i 

id 

i 

U 

-  X 

+ 

5 

cos 

0 

n 

sin 

e 

j 

j 

ji 

j 

ji 

j 

(2.14) 

V 

=  y 

+ 

K 

sin 

6 

+ 

n 

cos 

0 

i 

i 

ij 

i 

ij 

i 

V 

=  y 

+ 

5 

sin 

6 

+ 

n 

cos 

0 

j 

j 

ji 

j 

ji 

j 

Equation  2.13  restricts  the  motion  of  body  i  to  be  along  a  line  that  is 
at  a  specified  distance  from  and  parallel  to  the  line  of  relative  motion. 

The  second  scalar  equation  of  constraint  can  be  obtained  by  noting 
that  r  and  r  must  be  parallel.  This  is  true  since  both  of  these 

ij  ji  _ 

vectors  are  perpendicular  to  r  .  In  three  dimensions,  this  condition 

P 

can  be  expressed  as 

7  x  7  -0  (2.15) 

ij  ji 


Expanding,  the  component  perpendicular  to  the  x-y  plane  is 


♦  =(U  -x)(V  -  y  )  -  (V  -  y  )(U  -x)=0 

9  i  i  j  j  i  i  j  j 


(2.16) 
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Since  body  i  is  constrained  by  Eq.  2.13  to  move  parallel  to  r  , 

P 

Eq.  2.16  guarantees  that  body  j  will  also  move  along  a  line  parallel  to 

r  .  As  in  the  case  of  the  revolute  joint,  x  ,  y  ,  0  ,  x  ,  y  ,  and  0 
P  i  i  i  j  j  j 

are  the  state  variables  in  the  equations  of  constraint.  Variables  £  > 

ij 

n  ,  £  ,  and  n  and  related  to  the  design  variables. 

ij  ji  Ji 

2.2.3  Solution  Technique  for  State  Equations 
Constraint  equations  for  the  two  typical  joints  considered  here, 

Eqs.  2.10,  2.13,  and  2.16,  are  geometrically  nonlinear,  due  to  the  pres¬ 
ence  of  transcendental  functions  of  state  variables.  The  position 
state  equation  is  thus  nonlinear  and  a  solution  technique  that  is  appli¬ 
cable  to  nonlinear  equations  must  be  employed.  One  of  the  commonly  used 
techniques  for  solution  of  nonlinear  equations  is  Newton’s  method  [21]. 

Consider  the  position  state  equation  of  Eq .  2.6  for  the  entire 
system, 

$(z,b,a)  =  0 

Before  any  attempt  is  made  to  solve  this  nonlinear  system,  variables 
b  and  a  must  be  specified.  This  is  reasonable,  since  in  most  iterative 
design  algorithms  the  design  variable  b  is  estimated  before  the  synthesis 
procedure  is  initiated.  The  vector  of  input  parameters  a  is  generally  a 
part  of  the  problem  specifications.  The  only  unknowns  in  Eq.  2.6  are 
the  state  variables  z.  Equation  2.6  is  thus  a  system  of  3n  nonlinear 
equations  in  3n  variables  and  a  unique  solution  of  this  system  will 
exist  locally  if  the  stipulations  of  the  implicit  function  theorem  are 
satisfied;  ie. ,  if  the  matrix 
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9$(z,b,a) 

9z 


3<tl 


9z 


(2.17) 


| 3nx3n 

is  nonsingular.  This  condition  is  satisfied  if  the  system  of  constraints 
consists  of  no  redundant  joints. 

The  Newton  method  [21]  initially  requires  that  the  state  variable 
z  be  estimated.  The  method  then  computes  updates  Az  to  this  state  to 
obtain  an  improved  value  for  the  state.  The  improved  approximation  is 
given  by  [21] 


(i+1)  (D  (i) 

z  =  z  +  Az 

where  i  is  an  iteration  counter,  i  >  0, 


(2.18) 


Az 


(i) 


is  the  solution  of 


_i£  (z*,b,a)  Az*  =  -  $(z*,b,a)  (2.19) 

9z 

and  the  matrix  [9$/9z]  is  called  the  Jacobian  matrix  of  the  system. 

For  large  mechanisms,  the  system  of  Eq.  2.6  can  be  quite  large, 
giving  rise  to  a  large  system  of  linear  equations  in  Eq.  2.19.  Examin¬ 
ing  the  kinematic  constraint  equations  for  the  two  types  of  joints,  Eqs. 
2.10,  2.13  and  2.16,  it  can  be  concluded  that  these  pairs  of  constraint 
equations  involve  only  the  state  variables  of  the  two  bodies  that  they 
connect.  These  equations  are  thus  weakly  coupled  and  the  Jacobian 
matrix  in  the  left-hand  side  of  Eq.  2.19  is  highly  sparse.  Efficient 
sparse  matrix  codes  [22]  can  thus  be  used  for  solution  of  Eq.  2.19. 

Repeated  solution  of  systems  of  equations  similar  to  Eq.  2.19  are 
required  very  often  in  the  following  chapters.  To  perform  these 
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computations  efficiently,  the  sparse  matrix  code  initially  does  a 

symbolic  LU  factorization  of  the  coefficient  matrix.  Subsequent 

solutions  of  linear  systems  with  the  same  coefficient  matrix,  but  with 

different  right-hand  sides,  can  be  carried  out  very  efficiently. 

Since  mechanisms  are  being  synthesized  to  perform  over  a  range  of 

input  parameters  a,  solution  of  the  state  equations  would  be  required  at 

specific  values  of  a.  The  process  of  obtaining  the  solution  of  the 

position  state  equation  for  a  specified  value  of  a  can  be  repeated  to 

obtain  the  solution  for  a  desired  sequence  of  input  variables  a^ .  The 

numerical  efficiency  of  such  a  sequence  of  calculations  is  very  good 
j+1  j  j 

if  a  is  close  to  a  ,  since  z(a  )  serves  as  a  good  starting  estimate 

j+1  j+1  j 

in  the  computation  of  z(a  ).  If,  however,  a  and  a  are  not  close, 
then  an  update  6z(a^)  to  z(cf*)  is  required  to  produce  a  reasonable 
estimate  for  this  computation.  One  such  update  can  be  obtained  by 
linearizing  the  position  state  equation  of  Eq.  2.6,  keeping  b  fixed. 

9$(z»b>a3)  6z(aj)  +  9$(z>b>aJ)  =  0 

9z  9a 


9$(z,b,aJ)  6z  ( ^j )  .  _  3Kz>b,aj )  5a 

9z  9a 


(2.20) 


j+1  j 

where  6a  =  a  -  a  . 

j+1 

An  improved  estimate  for  z(a  )  can  thus  be  written  as 


0  j+1  j  j 

z  (a  )  =  z(a  )  +  6z(a  ) 


(2.21) 


where  6z(a  )  is  the  solution  of  Eq.  2.20.  Equation  2.20  has  the  same 


coefficient  matrix  as  Eq.  2.19,  so  its  numerical  solution  is  quite 


efficient. 


21 


2.3  Velocity  Analysis  Of  Mechanisms 
Most  mechanisms  are  driven  by  input  sources  that  give  input  links 
of  the  mechanism  finite  velocity.  It  is  then  necessary  to  determine 
velocities  of  the  remaining  links  in  the  mechanism,  during  motion  of  the 
mechanism,  over  the  prescribed  range  of  inputs. 

The  state  equation  for  the  mechanism,  Eq.  2.6,  is  required  to  hold 
for  all  time.  It  can,  therefore,  be  differentiated  with  respect  to  time 
to  obtain 

£_  $(z,b,a)  =  9$(z>b»ot)  z  +  3®(z»b,a)  „  =  0  (r) 

dt  3z  9^ - 


The  above  equation  can  be  rewritten  as 


_  3z  J 


z  =  -  3$(z,b,a)  a 
3a 


(2.23) 


3n 


where  z  e  r  is  the  vector  of  generalized  velocities  of  members  of  the 
system. 

Equation  2.23  is  the  velocity  state  equation.  This  equation  is 
linear  in  velocities  and  has  the  same  coefficient  matrix  as  Eq .  2.19. 

As  stated  in  section  2.2.3,  this  matrix  is  sparse  and  its  symbolically 
factorized  LU  form  has  been  determined  and  stored.  The  solution  of  Eq. 
2.23  is  thus  the  same  as  solving  Eq.  2.19,  with  a  different  right-hand 
side.  The  solution  of  this  equation  is  thus  very  efficient. 

The  components  of  the  vector  on  the  right  side  of  Eq.  2.23  are  not 
so  obvious.  Time  does  not  explicitly  appear  in  the  kinematic  constraint 
equations.  However,  in  the  driving  equations,  the  input  parameter  a 
could  be  given  as  an  explicit  function  of  time.  The  partial  derivative 
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of  Eq.  2.6  with  respect  to  time  can  be  written  as 


3$(z,b,a)  a  _ 
9a 


9$k(z,b)  * 

9a 

d  • 

9$  (z,b,a)  a 

3a 


(2.24) 


k 

Since  the  kinematic  constraints  $  -  0  do  not  depend  explicitly  on  a, 


Eq.  2.24  can  be  rewritten  as 


3S>(z,b,a)  a  _ 
3a 


0 

d 

3$  (z,b,a) 

3a  a 


(2.25) 


where  a  e  R  is  the  vector  of  first  time  derivative  of  input  parameters, 
d 

The  matrix  9$  (z,b,a)  in  Eq.  2.25  depends  on  the  form  of  the  driving 
9a 

constraints,  so  it  will  generally  be  problem  dependent. 


2.4  Acceleration  Analysis  of  Mechanisms 
Whenever  a  velocity  input  is  supplied  to  a  mechanism,  some  links 
experience  accelerations.  This  is  true  even  if  the  inputs  to  the 
mechanism  occur  at  constant  velocity.  Computation  of  accelerations  is 
important,  since  the  forces  experienced  by  links  in  the  mechanism  depend 
directly  on  acceleration. 

As  noted  in  Section  2.3,  the  position  state  equations  of  Eq.  2.6 
are  required  to  hold  over  the  entire  range  of  inputs,  so  the  velocity 
state  equation  of  Eq .  2.23  can  be  differentiated  once  again  with  respect 


to  time  to  obtain, 
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2 

—  *(z,b,a) 

dt2 


r  d$(z,b>a)~l  z  +  a 

p$(z,b,a)  1  5 

• 

i  +  ft  =  0 

L  9z  J  L3 

z  L  3z  J 

(2.26) 


Equation  2.26  can  be  written  as 


rM(z,b,ot)1  z  .  _  |"a  f^Cz^.g)  J  J  _  n  (2.27) 

L  3z  J  |J5z  L  3z  J  J 

where 

Q  =  J__  f~3$(z  ,b,  a)  '1  ‘a  +  3$(z,b>a)  0 
3z  L  3a  J  3a 


a  e  R  is  the  vector  of  second-time  derivatives  of  input  parameters 
3n 

and  z  e  R  is  the  vector  of  generalized  accelerations. 

Equation  2.27  is  the  acceleration  state  equation.  This  is  a  system 
of  linear  equations  with  the  same  coefficient  matrix  as  Eq.  2.19.  All 
the  desirable  properties  of  this  coefficient  matrix,  as  stated  in 
section  2.3,  still  hold.  The  solution  of  Eq .  2.27  is  thus  efficient. 

The  right-hand  side  of  Eq.  2.27  involves  z,  which  requires  that  the 
velocity  state  equation  of  Eq.  2.23  be  solved  before  Eq.  2.27  can  be 
solved.  The  second  term  in  the  right  side  of  Eq .  2.27  involves 
derivatives  of  constraints  with  respect  to  time.  Differentiating  Eq. 
2.25  with  respect  to  time  gives 


0 


- 

- 

3 

3$d(z,b,a)  . 

• 

34>  (z,b,a) 

•• 

3a 

a 

a  + 

a 

. 

3a  ^ 

3a 

_ 

(2.28) 
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CHAPTER  3 

FORCE  ANALYSIS  OF  MECHANISMS 

For  realistic  design  of  mechanisms,  it  is  necessary  to  impose 
stress  constraints  on  links  and  force  constraints  on  joint  bearings. 

This  requires  that  a  force  equation  be  derived  to  express  internal 
forces  on  links,  in  terms  of  the  externally  applied  forces  and  system 
velocities  and  accelerations.  The  applied  forces  could  be  forces  due 
to  gravity,  spring  damper  actuator  forces,  or  forces  from  other  external 
sources.  Two  approaches  are  taken  to  arrive  at  the  force  equations.  It 
is  also  shown  that  these  two  approaches  essentially  lead  to  the  same 
equations  for  equilibrium. 


3.1  Equilibrium  Equation  from  the  Principle 
of  Virtual  Work 

Figure  3.1  shows  a  body  with  a  body-fixed  coordinate  system  0  x  y  . 

An  externally  applied  force  F  and  an  external  moment  T  act  on  this 

ik  it 

body.  The  point  of  application  of  force  F  is  located  by  the  vector 
S  in  the  body-fixed  coordinate  system.  The  virtual  work  of  all 
external  forces  acting  on  body  i  can  be  written  as  [23] 


<5w  =  y  F  •  6(R  +  S  )  +  T  T  •  60  (3.1) 

i  k=l  ik  i  ik  zti  it  i 


where 


N  =  Total  number  of  forces  acting  on  body  i 
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•Figure  3.1  Force  Acting  on  Body  i 


26 


M  =  Total  number  of  moments  acting  on  body  i 
i 


Since  R  ,  S  ,  and  6  are  functions  of  the  position  state  variables  z, 
i  ik  i 

Eq.  3.1  can  be  written  as 


(3R  8S  \  Mi 

1  +  ik  Sz  +  l 

9z  9z  /  tr*l 


T 

i£ 


The  virtual  work  for  a  system  of  n  bodies  can  then  be  written  as 
n 

SW  ss  £  6w  (3.3) 

i=l  i 

Since  the  state  equations  for  position  of  Eq .  2.6  are  a  system  of 
workless  constraints,  the  principle  of  virtual  work  [23]  requires  that 
SW  of  Eq.  3.3  be  zero,  for  all  virtual  displacements  that  are  consistent 
with  the  constraints.  These  virtual  displacraents  are  all  Sz 
satisfying. 


6$ 


9$ 

9z" 


6z  -  0 


where  $  =  $(z,b,a). 


(3.4) 


3n 


Farkas  Lemma  [20]  now  guarantees  the  existence  of  a  vector  p  G  R 
of  multipliers  such  that 


T 

<5W  -  y  6$  =  0 


(3.5) 


3n 

for  all  <$z  e  R  .  Substituting  from  Eqs.  3.2,  3.3,  and  3.4  into 
Eq.  3.5,  one  has 


i 
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9R 

3S  \ 

£ 

90  ( 

1  + 

ik 

1+  l 

T  .  i 

■55“ 

-55-7 

1=  1 

ii  IzJ 

T 

y 


6z  =  0 


(3.6) 


which  must  hold  for  all  Sz  e  R  .  Since  6z  is  arbitrary  in  Eq .  3.6, 
each  of  its  components  can  be  varied  independently,  to  obtain  the  matrix 
equation 


I 


i=l 


F 

ik 


as 

ik 


T5~/ 


+ 


I1 

Jt=l 


T 

il 


Equation  3.7  can  be  rewritten  in  the  form 


(3.7) 


(3.8) 


where 


n 


H-  l 

i=l 


F 

ik 


/  9R  as  \ 
i  +  ^  + 

\'5i“  ~1Sz~  j 


I" 

1=1 


T 

il 


(3.9) 


Since  the  coefficient  matrix  of  Eq.  3.8  is  the  transpose  of  the  Jacobian 
matrix  of  Eq.  2.19,  this  linear  system  is  guaranteed  to  have  a  unique 
solution.  This  Jacobian  matrix,  as  stated  in  Section  2.2.3,  has  already 
been  symbolically  factored  and  stored.  Solving  Eq .  3.8  is  thus  the  same 
as  solving  Eq .  2.19,  but  with  a  transposed  coefficient  matrix  and  a 
different  right-hand  side.  The  solution  of  Eq.  3.8  is  thus  very 


efficient. 
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3.2  Force  Equations  From  Lagrange’s 
Equations  of  Motion 

Lagrange’s  equations  of  motion  for  a  dynamic  system  can  be  applied 
as  force  equations  for  kinematic  systems.  Considering  Lagrange’s 
equations  for  a  general  constrained  mechanical  system  [23], 


d 

dt 


jn_) 

/ 

KJ 

1 

2L  +  l  » 

8z,  *=1  * 


9z 


£  =  (Q  ) 
k 


+  (Q, ) 

nc  k  c 


k  1,...,  3n 


(3.10) 


where 


T 


(Q  ) 

k  nc 


=  kinetic  energy  of  the  system 

=  nonconservative  generalized  force  corresponding  to  the  kth 
generalized  coordinate 


(Q  )  =  conservative  generalized  force  corresponding  to  the  kth 

k  c 

generalized  coordinate 

r  =  total  number  of  constraints  on  the  system,  in  the  present 
context  r  =  3n 

y  =  vector  of  time-dependent  Lagrange  multipliers 

In  this  form,  Lagrange’s  equations  are  a  system  of  3n  equations 

in  3n  z  ’s  and  3n  y  ?s.  The  3n  equations  of  constraint  have  also  to  be 
k  i 

considered  along  with  Eq.  3.10,  to  solve  for  the  6n  unknowns  (3n  z  ’s  + 

k 

3n  y^’s).  However,  since  the  3n  equations  of  constraint  have  already 

been  solved  for  3n  z  ’s,  Eq.  3.10  has  only  3n  y  ’s  as  unknowns.  The 

k  i 
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generalized  velocities  and  accelerations  for  the  system  have  already 
been  determined.  Hence,  the  first  term  in  Eq.  3.10  is  a  known  quantity. 
Since  the  kinetic  energy  of  the  system  does  not  directly  depend  on  the 
generalized  coordinates,  the  second  term  in  Eq.  3.10  is  zero.  Equation 
3.10  can  thus  be  rewritten  as 
/9$  \ 

l  \  I — £)-«!>  +  (Q  )  -  k-1 . 3n 

*"  knc  kc  dt(3J 


The  right-hand  side  of  Eq.  3.11  may  be  written  as 
(Q  )  =  (Q  )  +  (Q  )  _  d  /9T 


nc 


k  c  dt 


(3.11) 


(3.12) 


9z 


k- 


where 


(Q  )  is  the  total  generalized  forces  corresponding  to  the  gener- 


ized  coordinate  z  . 

k 


Using  Eq.  3.12,  Eq.  3.11  can  be  written  as 

k  =  1,...,  3n 


3n 

l  V  _!  =  (Q  )  , 

4  9z-  k  s 


(3.13) 


*Zk 

Expressing  Eq .  3.13  as  a  matrix  equation,  gives 
T 

9$' 

—  I  y  =  G 
9z 


(3.14) 


where 


G  is  the  vector  of  (Q  )  ,  k  =  l,...,3n 

k  s 

The  force  equation  obtained  from  the  Lagrange  equations  of 
motion,  Eq.  3.14,  and  the  equilibrium  equation  obtained  from  the 
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principle  of  virtual  work,  Eq.  3.8,  are  of  the  same  form.  However,  the 
right  side  of  Eq.  3.14  includes  dynamic  effects.  These  dynamic  effects 
include  D'Alembert's  forces  and  moments  and  forces  and  torques  due  to 
spring  dampers.  For  mechanisms  being  acted  upon  by  static  forces  only, 
it  is  possible  to  show  that  the  right  side  of  Eq.  3.8  and  3.14  are 
equivalent. 

The  kinetic  energy  of  the  system  can  be  written  as 
3n  .  2 

T  =  £  m  (z  )  (3.15) 

k=l  k  k 

where 


m 

k 


=  mass  of  body  i  if 

k  =  {3i  -  2,  3i  -  1  |  1  <  i  <  n} 


=  mass  moment  of  inertia  of  body  i  if 


k  =  3i  ,  1  <  i  <  n 

Substituting  kinetic  energy  T  from  Eq.  3.15  into  the  third  term  of 
Eq.  3.11  gives 


d  /  3T  \  =  d  mz  =mz  (3.16) 

dt  (  T"  I  dt  k  k  k  k 


The  third  term  in  Eq.  3.11,  depending  on  index  k,  is  thus  D'Alemberts 

force  or  moment.  This  term  can  be  denoted  by  a  generalized  force  term 

(Q  )  .  The  right  side  of  Eq.  3.11  can  thus  be  written  as 
k  D 


(Q  )  -  {(Q  )  +  (Q  )  "  (Q  )  }  0*17) 

k  s  k  nc  k  c  k  D 

Consider  the  mechanical  system  being  acted  upon  by  a  system  of 

forces  ¥  and  moments  ¥  .  The  generalized  forces  due  to  this  force 

ik  ii 

system  can  be  directly  written  as  [24] 
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n  N.  9p  Mi  96 

(Q  )  =  l  l  F  •  i*+  [  T  •  i 

k  F  i-1  *- 1  U  1 r  A=1  H  ~te~ 

*—  k  k 


(3.18) 


where  P  is  the  vector  locating  the  point  of  application  of  force  F 

iH  a 

in  the  global  reference  frame. 

Since  the  generalized  coordinates  being  considered  in  Eq .  3.11  are 
the  same  as  those  in  Eq.  3.18,  the  corresponding  generalized  forces  can 
be  equated;  i.e. , 


(Q  )  -  (Q  )  -  (Q  ) 

kg  k  F  k  D 


(3.19) 


From  Fig.  3.1,  P  can  be  denoted  as 
i  i 


P  =  R  +  S 
i£  i  ±i 


(3.20) 


Substituting  for  the  right  side  of  Eq.  3.19  from  Eq.  3.18  gives 


n  %  /  9(R  +S  )\  Mi  39 

(q  )  =  y  y  f  .  I  *  **  1-+  y  ■  t  •  1  -  (q  )  (3. 21) 

k  s  i=l  Jt=l  3z  /  A=1  3z  k  D 

\  k  /  k 


Writing  the  above  equation  in  vector  form,  one  has 


n  N±  /  9(R  +S  )\  Mi_  36  _ 

G=y  l  F  £  T  •  i-(Q)  (3.22) 

i=l  £»1  i£  V  Sz  /  £=1  It  D 


By  directly  comparing  Eqs.  3.9  and  3.22,  it  can  be  seen  that  the 

right  sides  of  Eqs.  3.8  and  3.14  differ  by  the  term  (Q)  •  Since  (Q) 

D  D 

corresponds  to  the  D’Alembert  forces  or  moments,  for  static  mechanisms, 
Eqs.  3.8  and  3.14  are  completely  equivalent.  For  mechanisms  with 
dynamic  effects,  however,  only  Eq.  3.14  is  valid. 
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3.3  Relating  Lagrange  Multipliers  to  Joint  Reactions 

The  solution  of  Eq.  3,14  is  the  vector  of  Lagrange  multipliers  p. 
To  determine  joint  reactions,  and  the  subsequent  forces  in  the  members, 
it  is  necessary  to  relate  these  Lagrange  multipliers  to  the  joint  reac¬ 
tion  forces.  A  method  used  in  Ref.  13  will  be  used  here  to  develop 
relationships  between  the  Lagrange  multipliers  and  joint  reaction 
forces. 

Consider  a  mechanism  with  n  bodies  and  l  independent  joints,  with 
the  state  equations  for  position  given  by  Eq.  2.6.  For  a  typical  joint 
in  the  mechanism,  the  equations  of  constraint  are 

<f>  =0 

r 

(3.23) 

<f>  =0 

r+1 


In  differential  form,  Eq.  3.23  can  be  written  as 
3n  8<f> 

<5<j>  -  \  r  6z 

r  i=l  Jz~  1 

1  (3.24) 

3n  8(f) 

<S<J>  =  I  r+1  Sz 

r+1  i=l  dz  i 

1 


If  this  joint  were  to  be  ’broken' ,  the  equations  of  constraint  of 

Eq.  3.23  would  no  longer  hold.  However,  if  the  defects  (violations)  in 

constraint  equations  are  given  by  fid  and  fid  ,  respectively,  then  the 

r  r+1 

following  conditions  hold 


fi(J) 

+  fid 

=  0 

r 

r 

fi* 

+  fid 

=  0 

r+1 

r+1 

(3.25) 
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The  displacements  fid  and  fid  are  functions  of  the  generalized 

r  r+1 

coordinates  and  can  be  written,  s 
3n  3<|>  3n 

fid  =  -  l  r  §z  =  l  c  5z  (3.26) 

r  i=l  dz  i  i=l  1 

i 

3n  9<|)  3n 

fid  =  -  l  r+1  5Z  =  y  c  fiz  (3.27) 

r+1  i=i  9z  i  i=i  i,r+l  i 

i 

where 

c  and  c  are  functions  of  the  generalized  coordinates. 

i,r  i,r+l 

The  virtual  work  of  the  joint  reactions,  R  and  R  in  the  direc- 

r  r+1 

tion  of  d  and  d  ,  respectively,  can  be  written  as 
r  r+1 


6W  =  R  fid  +  R  6d 
k  r  r  r+1  r+1 


(3.28) 


where 

Joint  reactions  R  and  R  can  be  either  forces  or  moments,  and 
r  r+1 

k  is  the  number  of  the  joint  that  has  been  ’broken’. 

If  each  joint  and  driving  constraint  in  the  system  is  'broken' ,  the 
virtual  work  of  the  reaction  forces  is  given  by 


Z  2 Z+ m  3n 

fiW  =  y  5w  +  y  R  fid  =  y  R  fid 
R  k=l  R  p=2Z+l  P  P  r=l  r  r 


(3.29) 


where  p  is  the  index  of  the  driving  constraint.  Using  Eq .  3.25, 
Eq.  3.29  can  be  written  as 


3n 

5W  =  ~y  R  6(J> 

r  r_i  x  r 

Substituting  fromEq.  3.24, 

3n  3n  3  <j> 

5w  =  -  y  y  r  _l  Sz 

R  i=l  r=l  r  9z  i 

i 


(3.30) 


(3.31) 
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The  virtual  work  of  the  external  and  dynamic  forces  can  be  expressed 
in  the  form 


3n 

<SW  =  l  (Q  )  <Sz 
E  i=i  i  S  i 


(3.32) 


th 


where  (Q  )  is  the  generalized  force  associated  with  the  i  generalized 
i  S 

coordinate.  Since  the  joints  are  considered  broken,  the  equations  of 

constraint  of  Eq.  2.6  can  be  ignored.  Hence,  6z^  is  arbitrary  and 

independent.  Setting  the  total  virtual  work  to  zero,  6W  +  <$W  =0, 

R  E 

for  all  6z  gives 


5n  3<f> 

(Q  )  -  I  R  _£  , 

1  s  r=l  r  3z 

i 

In  matrix  form,  Eq.  3.33  becomes 


i  —  1,...,  3n 


(3.33) 


3$1TR=  G  (3.34) 

8z 

where 

R  =  vector  of  reaction  forces  R  and 

r 

G  =  vector  of  generalized  forces 

Comparing  Eqs.  3.34  and  3.14,  it  can  be  concluded  that 


R  =  u  (3.35) 

i  i 

Thus,  the  Lagrange  multipliers,  computed  as  a  solution  of  Eq .  3.14,  are 
the  reaction  forces  or  moments  in  the  joints  or  due  to  the  driving  con¬ 
straints.  The  definition  of  joint  reaction  forces  R  ,  however,  depends 
upon  the  definition  of  the  constraint  defects  6d^. 

Consider  now  the  two  typical  joints  described  in  Chapter  2.  The 
joint  reactions  in  the  revolute  joint  are  the  forces  in  the  global  X 
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and  Y  direction  respectively.  For  the  translational  joint,  the  joint 
reactions  are  the  moment  in  the  joint  and  a  force  acting  normal  to  the 
line  of  translation.  It  still  needs  to  be  shown  that  the  vector  of 
reactions  R  actually  represents  the  above-mentioned  reaction  forces. 


3.3.1  Reaction  Forces  in  Revolute  Joint 
Considering  a  revolute  joint  k  that  connects  body  i  and  body  j , 


Eqs.  2.7  can  be  rewritten  as 

(F  )  =[  R  +  F  -  R  -  r  | 

P  k  l  j  ji  i  ij/k 


(3.36) 


The  vector  r  can  be  written  in  terms  of  its  x  and  y  components, 
P 

r  and  r  as 
Px  Py 

(F  )  -  [r  r  ] 

p  k  P  n  i 

x  *y  k 

To  derive  the  equations  of  constraint  for  the  revolute  joint,  the 

condition  r  =  0  was  used  to  get 
P  ^  ^ 


(3.37) 


(r  )  = 

p  k  r 


In  differential  form,  when  <$r  *  0, 

P 


(3.38) 


Equation  3.38  implies  that  for  a  'broken'  revolute  joint,  the  variation 
of  the  x  and  y  constraint  equations  represents  the  negative  of  the  defect 


along  the  X  and  Y  axes,  respectively.  From  Eqs.  3.24  and  3.38, 
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6d  =  (fir  )  (3.39) 

r  Pxk 

6d  =  (fir  )  (3.40) 

r+i  Py  k 

where  index  r  corresponds  to  the  index  of  first  constraint  equation,  due 
to  revolute  joint  k. 

From  Eqs .  3.39  and  3.40,  it  can  be  concluded  that  the  Lagrange 
multipliers  corresponding  to  the  two  equations  of  constraint  of  a 
revolute  joint  are  the  x  and  y  reactions  in  the  revolute  joint. 


3.3.2  Reaction  Forces  in  Translational  Joint 
To  derive  the  reaction  forces  for  the  translational  joint,  a 
variation  of  the  procedure  derived  in  Section  3.3  is  used.  Since  the 
two  equations  of  constraint  for  any  joint  are  independent,  instead  of 
breaking  one  joint  at  a  time,  only  one  constraint  condition  for  the 
joint  is  broken  at  a  time.  The  unbroken  constraint  is  thus  still 
valid. 

Consider  the  second  equation  of  constraint  for  the  translational 
joint,  given  by  Eq.  2.15  as 

<J>  =  (U  -x  )(V  -y  )  -  (V  -y  )(U  -x  ) 

e  i  i  j  j  i  i  j  j 


In  differential  form,  this  is 

fid)  =  (5U  -fix  )(V  -y  )  +  (U  -x  )( 5V  -5y  ) 

6  i  i  j  j  i  i  j  j 

-  (6V  -fiy  )(U  -x  )  -  (V  -y  )(5U  -fix  )  (3.41) 

i  i  j  j  i  i  j  3 


Substituting  the  differential  form  of  U  ,  U  ,  V  ,  and  V  from  Eq.  2.14 

1  j  i  j 

into  Eq.  3.41  gives 
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6d>  =  {-(£  sin  6  +  n  cos  0  )(5  sin  0  +  n  cos  0  ) 

V9  1  ij  i  ij  i  ji  J  n  3 


-(£  cos  0  -  n  sin  0  )(5  cos  0  -  n  sin  0  ) } 

ij  i  ij  i  ji  J  Ji  3 


x  (  60  -  60  )  (3.42) 

i  3 


Equation  3.42  can  be  written  as 

6 d>  =  c  (69  -  60  )  (3.43) 

9  eu  i  J 


where  c  represents  the  term  in  curly  brackets  in  Eq.  3.42. 

0 

ij 

Thus  the  condition  <|>  =  0  constrains  relative  rotation  of  bodies  i  and  j. 

0 

With 

6d>  +  6d  =  0 
0  0 


as  in  Eq.  3.25 

6d  =  -  c 
0 


> 


(  60  -  60  ) 
i  j 


which  is  proportional  to  a  virtual  relative  rotation  of  bodies  i  and  j. 


This  shows  that 


p  6d  =  -p  c  (60  -  60  ) 

0  0  0  0.  .  i  i 

iJ  J 


SO  p  c 
0  0 


ij 


is  the  reaction  torque  acting  between  bodies  i  and  j,  due  to  a 


translational  joint. 


Consider  now  the  first  equation  of  constraint  for  the  translational 


joint,  given  by  Eq.  2.13  as 

<j>  =  (U  -x  )(U  -U  )  +  (V  -y  )  (V  -V  ) 

n  iiij  i  i  i  j 


In  differential  form,  this  is 

<5<j>  =  (<5U  -6x  )(U  -U  )  +  (U  -x  )(6U  -6U  ) 

n  iiij  ii  ij 


-  (6V  -6y  )(V  -V  )  +  (V  -y  )(6V  -6V  )  (3.45) 

iiij  ii  ij 


Substituting  the  differential  forms  for  U  ,  V  ,  V  and  V  from  Eq.  2.14 

i  j  i  j 

into  Eq.  3.45  gives 

<$<(>  =  [~(V  -y  )  66  ](x  -x  )  +  [(U  -x  )<S6  ](y  -y  ) 

n  iiiij  iij'ij 

+  [ (U  -x  )  ](6x  -6x  )  +  [v  -y  ](5y  -5y  ) 
i  i  i  j  iiij 


-  (V  -y  )(-£  cos  0  +  q  sin  0  )60 

i  i  ji  j  ji  j  i 


-  (5  cos  0  -  q  sin  0  )(£  sin  0  +  q  cos  0  )60 

ij  i  ij  i  ij  i  ij  i  i 


+  (U  -x  )(£  sin  0  +  q  cos  0  ) 60 

i  i  ij  i  ij  i  j 


“  (5  sin  0  +  q  cos  0  )(£  cos  0  -  q  sin  0  ) 60 

ij  i  ij  i  ji  j  ji  j  j 


The  above  expression  can  further  be  simplified  to 

6 4>  =  [-(V  -y  )  60  ](x  -x  )  +  [(U  -x  )  60  ](y  -y  ) 

n  iiiij  iijij 

+  [u  -x  ](6x  -6x  )  +  [v  -y  ] (  6y  —  6y  ) 
i  i  i  j  iiij 

+  {(5  sin  0  +  q  cos  0  )(£  cos  0  -  q  sin  0  ) 

ij  i  ij  i  ji  j  ji  j 

-  (5  cos  0  -  q  sin  0  )(£  sin  0  +  q  cos  0  ) } (  60  -60  ) 

ij  i  ij  i  ij  i  ij  i  '  i  j 
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Since  each  of  the  two  equations  of  constraint  for  this  joint  are  being 
'broken'  one  at  a  time,  the  second  equation  of  constraint  for  this  joint, 
Eq.  2.16,  is  valid  here.  From  Eq.  3.43,  this  implies 

6<J>  =  c  (60  -  60  )  =  0  (3.47) 

9  9u  1  j 


Coefficient  c  has  to  be  nonzero  to  get  a  nontrivial  constraint 

9ij 

equation.  Equation  3.27  thus  implies 

(60  -  66  )  =  0  (3.48) 

i  j 

Substituting  Eq.  3.48  into  Eq.  3.46  gives 

6i{)  =  [-(V  -y  )60  ](x  -x  )  +  [(U  -x  )60  ](y  -y  ) 

n  iii  ij  iijij 

+  [u  -x  ](6x  -6x  )  +  [v  -y  ](6y  -6y  )  (3.49) 

i  i  i  j  i  i  i  j 


To  develop  an  interpretation  of  Eq.  3.49,  consider  Fig.  2.3 


again.  Write  a  unit  vector  n  along  r  as 
1 

n  =  { (U  -x  )  T  +  (V  -y  )  J} 

If  I  1  i  i  i  i 

ij 


(3.50) 


Since  d>  is  still  a  valid  constraint ,  n  is  also  parallel  to  r  .  The 

0  ji 

projection  of  the  position  vectors  of  the  body-fixed  coordinate  systems 


of  bodies  i  and  j  along  n  can  be  expressed  by  n^  and  n  ,  respectively  as. 


n  =  n  •  R 
i  i 


n  =  n  •  R 

j  j 


(3.51) 

(3.52) 


Denote  n  as 


n  =  n  I  +  n  J 
1  2 


(3.53) 
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where 

n  =  (U  -x  ) / | r  | 

1  i  i  ij 

n  =  (V  -y  )/|r  | 

2  i  i  ij 

In  differential  form,  Eqs.  3.51  and  3.52  can  be  written  as 


fin 

-  fin  x  +  n  fix  + 

6n  y 

+  n  fiy 

(3.54) 

i 

1  i  1  i 

2  i 

2  i 

fin 

=  fin  x  +  n  fix  + 

6n  y 

+  n  fiy 

(3.55) 

j 

i  j  i  j 

2  j 

2  j 

The  virtual  displacement  of  body  i  relative  to  body  j  along  the 


direction  of  vector  n  can  be  expressed  as 


Sn  -  6n  =  6n  (x  -x  )  +  6n  (y  -y  ) 
i  j  1  i  j  2  i  j 


+  n  ( 6x  -  8x  )  +  n  (  6y  -6y  )  (3.56) 

1  i  j  2  i  j 


Observing  that  |r  I  is  invariant  under  a  rotation,  the  differential 

ij 

forms  of  n  and  n  can  be  written  as 
1  2 


(y  -y  ) 

«Sn  =  -  i  i  66  (3.57) 

1  i 
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6n  -  6n  =  1  {[“(V  -y  ) 66  ](x  -x  ) 

i  i  j 


+  [u  -x  ]( 6x  -  6x  ) 

L  i  iJ  i  j 

(3.59) 


Comparing  Eqs.  3.49  and  3.59  gives 

6<j>  =  |  r  |  (  Sn  -  Sn  )  (3.60) 

n  ij  i  j 


Rewrite  Eq.  3.60  as 


6 c(>  =  c  (Sn  -  Sn  )  (3.61) 

n  n  i  j 

ij 

where  c  «  |7  | .  Thus  the  condition  $  =  0  constraint  relative 

ni:j  ij  n 

translation  of  bodies  i  and  j  in  the  direction  n. 

With 

6<j>  +  Sd  =0 
n  n 


as  in  Eq .  3.25, 

6d  =  -  c  (Sn  -Sn  ) 

n  "ij  1  j 

is  proportional  to  a  normal  virtual  relative  displacement  of  bodies  i 


+  [(iw^eJCypO 
+  )} 


and  j.  This  shows  that 


—  y  Sd  =  -y  c  (  Sn  -  Sn  ) 


n  n 


n  n, 


so  y  c 
n  n 


ij  1  J 

is  the  normal  reaction  force  acting  between  bodies  i  and  j, 


ij 


due  to  a  translational  joint 
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CHAPTER  4 

OPTIMAL  DESIGN  OF  MECHANISMS 

4.1  Introduction 

/ 

Chapters  2  and  3  provide  the  theory  necessary  to  compute  the 
kinematics  of  the  mechanism  and  the  forces  acting  on  the  links  of  the 
mechanism.  It  should,  therefore,  be  possible  to  put  design  constraints 
(bounds)  directly  on  these  variables  or  on  functions  of  these  variables. 
Also,  it  should  be  possible  to  make  extreme  any  of  the  state  variables 
or  their  functions,  subject  to  constraints. 

Most  extreme  algorithms  used  require  that  explicit  derivatives, 
with  respect  to  design  variables,  of  the  cost  and  constraint  functions 
be  provided.  Computing  derivatives  of  functions  involving  only  design 
variables  is  easy.  However,  for  functions  involving  state  variables, 
the  dependence  on  design  arises  indirectly,  through  the  state  equation. 
Derivative  computation  of  such  functions  is  thus  not  as  simple  as  that 
for  explicit  functions  of  design  variables.  An  adjoint  variable  tech¬ 
nique  that  has  been  used  for  design  sensitivity  analysis  of  structural 
system  [20],  is  used  in  this  chapter  for  design  sensitivity  analysis  of 
mechanisms . 


4.2  Statement  of  the  Optimal  Design  Problem 


4.2.1  Statement  of  Continuous  Optimization  Problem 

A  general  class  of  optimal  design  problems  can  be  stated  as: 

S 

Find  a  design  b  R  to  minimize  the  cost  function 
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(4.1) 

(4.2) 

(4.3) 

(4.5) 

(4.6) 


(z,z,z,y,b)  =  0  ,  i  =  p  +  1»***>P  +  Q  (^*7) 

i 

Equations  4.1  to  4.7  define  a  general  optimal  design  problem.  Any 
type  of  design  constraint  can  be  treated  in  this  formulation,  as  long  as 
it  can  be  put  in  a  form  such  as  Eqs.  4.6  or  4.7.  The  representation  of 
the  cost  function  in  the  form  of  Eq.  4.1  does  not  restrict  the  technique 
from  being  applied  to  cost  functions  involving  state  variables.  An 

upper  bound  technique  [20]  that  is  used  in  such  cases  is  now  illustrated 

A  • 

for  a  general  function  (z,z,z,u,b),  the  maximum  value  of  which  is 
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required  to  be  minimized,  over  a  specified  range  of  input  parameters  a. 
Thus, 

A  A  ^ 

min  i|>  =  min  {max  ^  (z,z,z,y,b)}  (4.8) 

b  0  b  a  0 

The  above  formulation  of  the  cost  function  is  natural  for  kinematic 

optimization.  Since  state  variables  take  on  different  values  over  the 

entire  range  of  input  parameters,  it  is  natural  to  minimize  the  maximum 

value  of  functions  of  these  variables. 

Equation  4.8  represents  a  mini -max  problem  [20]  and  is  not  simple 

to  deal  with  directly.  A  scheme  commonly  employed  is  to  introduce  an 

A  •  “ 

artificial  design  variable  b  to  be  an  upper  bound  of  ^  (z,z,z,u,b). 
Therefore,  the  minimization  problem  in  Eq,  4*8  can  be  written  as 

■dn  ♦„ s  Vi  (4-9) 

subject  to  the  additional  constraints 


•  • 

(i)  (z,z,z,y,b)  =  {max  (z,z,z,y,b) }  -  b  <  0 


(4.10) 


(ii)  State  Equations  and  other  design  constraints. 

The  minimization  problem,  as  stated  in  Eqs.  4.9  and  4.10,  amounts  to 

A 

generating  a  minimizing  sequence  of  upper  bounds  of  the  function  tj>  . 

0 

Composite  design  constraints  in  the  form  of  Eqs.  4.6  and  4.7  are 


required  to  hold  over  the  entire  range  of  specified  input  parameters. 

Such  constraints  are  called  parametric  constraints.  Techniques  for 
making  cost  functions  extreme  which  are  subject  to  parametric  constraints 
have  been  developed  earlier  [25] .  Since  these  techniques  require  a  con¬ 
siderable  amount  of  additional  computation,  a  simpler  approximate  tech¬ 
nique  is  resorted  to.  The  range  of  input  parameters  a  is  discretized 
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into  a  set  of  grid  points  a  ,  j  =  1,...,t.  The  composite  design 
constraints  are  then  required  to  hold  at  every  point  on  this  grid. 


4.2.2  Statement  of  Discretized  Optimal 
Design  Problem 

The  optimal  design  problem,  with  a  grid  imposed  on  the  range  of  the 

input  parameters,  can  be  stated  as  follows: 

s 

Find  a  design  b  R  to  minimize  the  cost  function 

ip  =  if>  (b)  (4.11) 

0  0 


subject  to  State  Equations; 

(A.l)  Position  State  Equations  of  Eq.  2.6 

$(z^,b,a^)  =  0  ,  j  =  1 ,  • .  • ,  T  (4.12) 

where  t  =  number  of  grid  points  on  the  range  of  the  input  parameters. 
(A. 2)  Velocity  State  Equation  of  Eq.  2.23 


3$(zj ,b,a^  ) 


9z 


J  ,  j-1 . T  <4->3> 

8a 


(A. 3)  Acceleration  State  Equation  of  Eq.  2.27 


9$(z^  ,b,a*^ ) 

8 

z  =  - 

8$(z^  ,b,a^)  *j 

9z 

8z  J 

(A. 4)  Force  Equation  of  Eq.  3.14 

T 


8$(z^ ,b,a^ ) 


8z 


^  =  cP  > 


and  Composite  Design  Constraints; 
(B.l)  Inequality  Constraints 


z^  -  »  j  =  1 , . . . ,  t 

(4.14) 


j  -  1,...,T 


(4.15) 


j  *j  **j  j 

(z  ,z  ,z  ,ti  ,b)  <0,  i  =  l,...,p,  j  =  1,...,t  (4.16) 

i 
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(B.2)  Equality  Constraints 


,  j  *3  J  3  ,  . 
ip (z  ,Z  ,z  ,\1  ,b)  =  0, 


i  =  P+1, 


jP+H  »  3  _  1 » • « • » t 

(4.17) 


4.3  Design  Sensitivity  Analysis 
As  noted  in  the  introduction  to  this  chapter,  all  design 
optimization  algorithms  require  that  derivatives  of  cost  and  constraint 
functions,  with  respect  to  design  variables,  be  provided.  Having  stated 
the  optimal  design  problem,  it  is  now  possible  to  proceed  to  deriving 
the  derivatives  of  the  cost  and  constraint  functions.  This  derivation 
is  restricted  to  the  discretized  optimal  design  problem. 

The  first  variation  of  the  cost  function  of  Eq .  4.11  can  be  written 
simply  as 

3ip  (b)  T 

=  0  6b  =  iP  6b  (4.18) 

0  3b 


Since  the  cost  function  can  always  be  reduced  to  a  function  of  design 
variables  by  the  method  explained  in  Section  4.2.1,  the  sensitivity  of 
the  cost  function  can  very  simply  be  written  in  the  form  of  Eq.  4.18. 

The  first  variation  of  the  composite  design  constraint  can  be 
written  as 


(4.19) 
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Moving  all  terms 


in  Eqs.  4.22  through  4.25  involving  fib  to  the 


(4.23) 


(4.24) 


(4.25) 


right 


side,  this  system  of  equations  can  be  written  in  matrix  form  as 
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(4.26) 
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where 

j  ,  j  *J  “j  J  x 

A  =  A(z  ,z  ,z  ,p  ,b) 

J  .  j  *J  “j  j  v. 

B  =  B(z  ,z  ,z  ,y  ,b) 

j 

In  deducing  Eq .  4*28  from  Eq .  4*27,  it  is  required  that  the  matrix  A  be 
invertible.  The  existence  of  such  an  inverse  is  proved  in  Appendix  A. 
Equation  4.28  expresses  the  variation  of  state  variables  in  terms 

j 

of  variations  in  the  design  variables.  Substituting  SU  from  Eq.  4.28 
into  Eq.  4.20  gives 


(4.29) 


To  avoid  calculation  of  (A  )  , 

j  "I 

constraint  derivatives  and  (A  )  is 

.  .T 
ij 

vector  A  ,  i.e.  , 


the  product  of  the  row  vector  of 
denoted  by  a  composite  adjoint 


8  ij; 

a*. 

3^ 

i 

i 

1 

i 

dz 

• 

- 

3y 

3z 

9z 

(4.30) 


Equation  4.30  can  be  rewritten  as 


(4.31) 


The  composite  adjoint  vector  A  , 
12n  x  1  vector.  This  vector  can 

ij 

four  (3nxl)  adjoint  vectors  A^ 


as  represented  in  Eq.  4.30,  is  a 

be  written  as  a  composite  vector  of 
ij 

to  A  .  Equation  4.31  can  thus  be 


expanded  as 
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Equation  4.32  is  a  system  of  12n  linear  equations.  Rather  than  solve 
this  as  a  coupled  system  of  12n  equations,  the  form  of  the  coefficient 
matrix  makes  it  possible  to  solve  four  separate  linear  systems  of  3n 
equations  each. 

As  can  be  seen  from  Eq.  4.32,  there  is  no  coupling  between  the 
adjoint  vector  X  ^  and  any  of  the  other  adjoint  vectors.  Thus  the  last 


3n  equations  can  be  written  as 


which  has  the  same  coefficient  matrix 
2.2.3,  the  solution  of  linear  systems 
very  efficient,  since  the  LU  factored 
already  been  computed  and  stored. 

Equations  (6n  +  1)  to  (9n)  of  Eq 


r  i 

T 

ij 

/  \T 
-  /  3G\ 

ij 

~9i|; 

3$ 

X 

X 

i 

3z 

3 

~ 

4 

•• 

j 

\3z/j 

3z 

Since  X  has  already  been  determined 
4 

Eq.  4.34  can  be  rewritten  as 


3$ 

T 

Xij  - 

M’ 

3z 

3 

j 

i 9z  Jj 

1 4 

The  coefficient  matrix  in  Eq .  4.35  is 
matrix  of  Eq.  2.19.  The  remarks  made 
efficiency  of  Eq .  4.33  are  also  valid 


(4.33) 

as  Eq .  2.19.  As  noted  in  Section 
with  this  coefficient  matrix  is 
form  of  the  coefficient  matrix  has 

4.32  can  then  be  written  as 

(4.34) 

as  the  solution  of  Eq.  4.33, 

(4.35) 

the  transpose  of  the  coefficient 
above  about  the  solution 
for  Eq .  4. 35. 
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Continuing  this  process  of  backward  solution,  equations  (3n+l)  to 
(6n)  of  the  linear  system  of  Eq.  4.32  can  now  be  written  as 


+  Xij  + 

3 


X^ 

4 


(4.36) 

ij  ij 

where  X  and  X  are  the  solutions  of  Eqs.  4.35  and  4.33, 

3  4 

respectively . 

As  in  the  case  of  Eq .  4.35,  the  coefficient  matrix  of  Eq .  4.36  is 
the  transpose  of  the  coefficient  matrix  of  Eq.  2.19,  so  the  remarks  made 
about  the  solution  efficiency  of  Eq .  4.33  are  also  valid  for  Eq.  4.36. 

The  first  3n  equations  of  the  linear  system  in  Eq .  4.32  can  now  be 
written  as 
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where  X  ,  X  ,  X  are  the  solutions  of  Eqs.  4.36,  4.35  and  4.33, 

2  3  4 

respectively.  The  coefficient  matrix  of  Eq.  4.37  is  the  transpose  of 
the  coefficient  matrix  of  Eq.  2.19,  so  the  remarks  made  about  the 
solution  efficiency  of  equations  with  this  coefficient  matrix  are  also 
valid  for  Eq.  4.37. 

Equations  4.33  and  4.35  to  4.37  are  four  adjoint  equations  associ¬ 
ated  with  the  four  state  equations.  Now  that  the  solution  of  the 

•  .T 
ij 

composite  adjoint  vector  A  is  known,  Eq.  4.30  can  be  substituted 
into  Eq.  4.29  to  give 


(4.38) 


Equation  4.38  expresses  the  variation  of  a  composite  design  con¬ 
straint  in  terms  of  the  variation  in  design  only.  Comparing  Eqs.  4.21 

and  4.38,  it  can  be  concluded  that  the  quantity  premultiplying  6b  in 

th  < 

Eq.  4.38  is  the  design  sensitivity  vector  of  the  i  constraint,  at  the 
th 

j  grid  point.  Therefore, 


T  T 

*ij  =  tf-3  jjJ 


Substituting  for  B 


from  Eq.  4.26,  this  can  be  explicitly  written  as 
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Equation  4,39  gives  the  design  sensitivity  of  a  composite  design 
constraint  i  at  grid  point  j,  in  terms  of  derivatives  of  the  position 
state  equation  and  the  adjoint  variables. 

4.4  Design  Optimization  Algorithm 
4.4.1  Active  Set  Strategy 

With  the  design  sensitivity  information  computed  in  the  previous 
section,  one  can  proceed  to  implement  the  optimization  algorithm  of  his 
choice.  The  gradient  projection  algorithm  with  constraint  error 
correction  [20]  has  been  used  in  the  past  in  related  applications 
[26,  27].  This  algorithm  is  used  here  for  design  optimization.  Before 
implementing  this  algorithm  for  the  present  application,  some  observa¬ 
tions  can  be  made  that  will  enhance  efficiency  of  computation. 

In  connection  with  the  composite  design  constraints  of  Eqs.  4.16 
and  4.17,  it  can  be  observed  that  many  design  constraints  can  be  put  in 
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a  simpler  form  than  that  indicated  in  these  equations*  Constraints  that 
do  not  involve  the  state  variables  do  not  explicitly  or  implicitly  de¬ 
pend  on  the  input  variables  a*  Such  constraints  are  called  ’non¬ 
parametric 1  constraints*  Design  constraints  that  involve  the  state  var¬ 
iables  are  called  ’parametric1  constraints*  It  is  necessary  to  make 
this  distinction,  since  only  the  parametric  constraints  are  required  to 
be  satisfied  over  the  entire  range  of  input  variables  a*  For  a  given 
design,  the  nonparametric  constraints  need  only  be  evaluated  once. 
However,  parametric  constraints  must  be  evaluated  at  all  points  on 
the  interval  of  the  input  variables. 

An  active  set  strategy  may  be  adopted  to  determine  the  reduced  set 
of  active  constraints.  Since  equality  constraints,  parametric  or 
nonparametric,  are  always  active,  they  are  always  included  in  the 
active  set  Nonparametric  inequality  constraints  that  are  e-active 
are  also  included  in  the  reduced  set  i|>.  Parametric  inequality 
constraints,  due  to  their  dependence  on  the  input  variable  a,  must  be 
evaluated  on  all  points  on  the  grid  of  input  variable  a.  Some 
computational  efficiency  can  be  realized  from  the  fact  that  the  gradient 
projection  algorithm,  as  stated  in  the  following  section,  allows  only 
small  changes  in  design  leading  to  small  changes  in  state.  A  design 
constraint  with  a  large  violation  at  a  given  design  iteration  may  not  be 
fully  satisfied  at  the  subsequent  design  iteration.  This  is  so  because 
the  algorithm  used  for  optimization  here  uses  only  first-order 
information  about  the  design  constraints.  However,  design  constraints 
are  generally  nonlinear.  The  regions  of  the  input  variable  grid  in 
which  a  parametric  inequality  constraint  is  active  are  also  not  expected 
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to  change  rapidly  from  design  iteration  to  iteration.  It  is  thus 
possible  to  avoid  evaluating,  for  a  few  design  iterations,  a  parametric 
constraint  in  the  region  in  which  it  is  not  E-active. 

Design  iterations  in  which  the  parametric  constraint  is  evaluated 
at  each  point  on  the  a  grid  are  defined  as  ’sweep’  iterations. 

Iterations  in  which  the  parametric  constraint  is  evaluated  only  in  the 
active  region  are  called  ’nonsweep’  iterations.  The  interval  between 
two  ’sweep’  iterations  depends  upon  how  rapidly  the  active  regions  are 
changing.  For  ’nonsweep’  design  iterations,  it  is  not  necessary  to 
solve  the  state  equations  on  the  entire  range  of  the  input  variable. 
Considerable  computational  saving  can  be  realized  by  having  a  large 
number  of  ’nonsweep’  iterations  between  successive  ’sweep’  iterations. 
Since  new  active  regions  can  only  be  detected  during  ’sweep’  iterations, 
having  a  large  number  of  ’nonsweep1  iterations  separating  two  ’sweep’ 
iterations  could  lead  to  new  active  regions  going  undetected  for  a 
number  of  design  iterations. 

Two  alternative  definitions  could  be  used  to  define  active  regions 

on  the  grid  of  the  input  variable.  The  first  definition,  as  shown  in 

Fig.  4.1,  involves  determining  the  e-active  relative  maxima  of  the 

constraint  function  on  the  a-grid.  The  active  region  is  then  defined  to 

be  the  set  of  points  at  which  the  relative  maxima  occur  and  one  grid 

point  on  either  side  of  these  grid  points.  The  active  region  can  thus 

be  defined  by  the  index  set  I  , 

R 

I  «<I  U  I  U  I 
R  K  *M  \ 


(4.40) 


Constraints 
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where 


I  =  {j-1,  j,  j+1  |  *1J  >  -  e,  i  =  2  <  j  <  t-1; 

is  a  relative  maximum  point  } 

r  i 

I  -  j.  j+1  *  >  -  e,  i  =  j  =  1; 

*L 

is  a  relative  maximum  point  } 

r  i 

I  =  {j-1,  j,  I  ♦  >  “  e,  i  =  j  =  t; 

is  a  relative  maximum  point  } 

The  second  definition  of  active  region  involves  determining  all  the 

points  on  the  a  grid  at  which  the  constraint  function  is  e-active.  This 

set  is  defined  to  be  the  e-active  region.  Denoting  this  set  as  I  , 

E 

I  =  {j  |  ty2  >  -  e,  i  =  l,...,p,  1  <  j  <  t}  (4.41) 

E 

The  relative  maximum  strategy  has  been  suggested  for  defining 
active  regions  for  parametric  constraints  in  Ref.  [25].  When  applied  to 
constraints  in  mechanism  optimization,  this  strategy  many  times  causes 
rapid  oscillation  of  the  relative  maximum  point  on  the  a-grid.  A  switch 
to  the  e-active  strategy  overcomes  this  problem.  As  is  evident  from 
Figs.  4.1  and  4.2,  there  is  a  penalty  to  be  paid  for  this  switch,  since 
the  latter  strategy  requires  a  larger  number  of  grid  points  to  be 
included  in  the  e-active  region. 
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4.4.2  Gradient  Projection  Algorithm 

The  Gradient  Projection  Algorithm  [20]  for  design  optimization  can 

now  be  stated  in  the  following  steps : , 

0 

Step  :  Estimate  a  design  b  ,  and  impose  a  grid  on  the  range  of 
the  input  parameters . 

Step  2:  Solve  the  state  equations  of  Eqs.  4.12  to  4.15  for  state 
j  .j  ”1  3 

variables  z  ,  z  ,  z  ,  y  ,  respectively,  where  j  =  1,...,t  if  the 

current  iteration  is  a  'sweep'  design  iteration  or  j  I  or 

R 

I  .  Note,  the  first  design  iteration  must  be  a  sweep  iteration. 

G 

Step  3:  Determine  the  active  region,  depending  on  the  strategy 
chosen,  and  form  a  reduced  vector  i|>  consisting  of  all  e-active  non- 
parametric  inequality,  and  all  equality  parametric  and  nonpara- 
metric  constraints.  For  inequality  parametric  constraints,  the 


constraints  evaluated  in  the  active  regions  are  included  in  i|>. 

Step  4:  Compute  adjoint  variables  X  X  X  \  and  X  ^  from  Eqs. 

4  3  2  1 

4.33,  4.35,  4.36  and  4.37,  respectively,  and  construct  design 

ij 

sensitivity  vectors  Z  of  Eq.  4.39  for  the  constraint  functions  in 

ij  ij 

•  Form  the  matrix  A=[Jl  ]  ,  whose  columns  are  the  vectors  Z 


corresponding  to  constraint  functions  in  ip.  Thus,  6 xp  =  A  6b 


Step  _5:  Compute  the  vector  M  and  matrix  M  from  the  following 
relations; 

~T  -1  0 

M  =  A  W  Z  (4.42) 

0 

where  Z  is  defined  in  Eq.  4.18 

~T  -1  - 

M  =  A  W  A 

W 


(4.43) 
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Step  6:  In  the  first  iteration,  compute  a  parameter  y  related  to 


step  size,  as 


0T  -1  0 

£  W  i  (4.44) 

Y  =  0 

28  ip  (b  ) 

0 


where  8  is  the  desired  fractional  reduction  in  the  cost  function. 


The  usual  range  of  8  is  0.03  to  0.15.  In  succeeding  iterations. 


the  factor  y  is  adjusted  to  enhance  convergence  of  the  algorithm. 
“1  ~2 

Step  7:  Compute  p  and  p  from  eqs.  4.45  and  4.46, 


“1 

M  p  =  -  M  (4.45) 

M  Mo 


~2  ~ 

M  p  =  -  A\p  (4.46) 

n 


where  AiJ>  =  C  ip  ,  C  is  the  fraction  correction  of  the  constraint 

'P  ip 

desired,  usually  in  the  range  0.30  to  1.0. 

1  2 

Step  8:  Compute  6b  and  6b  from  Eqs.  4.47  and  4.48  as 


1  -1  fo  ~f| 

6b  =  W  U+ApJ  (4.47) 

2  -1  “2 

6b  =  -W  A  p  (4.48) 


Step  9:  Compute  an  update  in  design  6b  from  Eq .  4.49 

1  2 

6b  =  -  _J_  6b  +  6b  (4.49) 

2Y 

Step  10:  Update  the  estimate  of  the  optimal  design  using  Eq.  4.50 


1  0 

b  =  b  +  6b 


(4.50) 
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Step  11:  If  all  constraints  are  satisfied  to  within  the  prescribed 
tolerance  and 


1 16b1!  I  - 


S  1  2 

I  W  (6b  ) 
i=l  i  i 


1/2 


<  6 


(4.51) 


terminate  the  process.  Where  6  is  a  specified  convergence 
tolerance.  If  Eq .  4.51  is  not  satisfied,  return  to  Step  2. 
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CHAPTER  5 

NUMERICAL  EXAMPLES 

5.1  Example  1  -  Kinematic  Synthesis 
of  a  Path  Generator 


5.1.1  Problem  Description 

A  segment  of  a  straight  line  is  required  to  be  generated  by  a  point 

P  on  the  coupler  of  the  four-bar  mechanism  shown  in  Fig.  5.1.  In 

addition  to  having  the  lengths  of  various  links  as  design  variables 

(b^  to  b^),  the  orientation  of  the  base  link,  body  1,  is  a  design 

variable  (b^_).  The  orientation  of  the  reference  line  with  respect  to 

the  base  link,  about  which  the  input  variable  a  is  measured,  is  also  a 

design  variable  (b  ) .  The  other  two  design  variables  (b  and  b  )  are  as 
4  7  8 

indicated  in  Fig.  5.1.  The  length  of  the  base  link  is  kept  fixed  at  10 
units.  This  problem  is  the  same  as  Example  3  in  Ref.  26. 

5.1.2  Problem  Formulation 

The  motion  specification  for  this  problem  requires  that  the  devi¬ 
ation  of  the  y  coordinate  of  coupler  point  P  from  zero  in  the  global 
coordinate  system  be  minimized.  The  position  vector  for  point  P,  in 
terms  of  design  and  state  variables,  can  be  written  as 

b2 

Rp  =  x3  +  cos  0g( - 2  +  by  cos  bg)  -  bySin  bgSin  6^  I 

b? 

=  y3  +  sin  6g( - J  +  b7  cos  b8^  +  b7cos  03sin  bg  J  (5.1) 
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Equation  5.1  can  be  symbolically  expressed  as 


R  =  (R  )  I  +  (R  )  J 


(5.2) 


where  Rp  and  Rp  represent  the  x  and  y  coordinates  of  point  p  in  the 
Fx  *y 

global  coordinate  system. 


An  artificial  design  variable  b  is  introduced  such  that 

9 


<  bQ  ,  a 


„  .  <  a  <  a 

min  max 


where  the  range  of  a  is  given  as  a  =  -0.3491  radians  and  »  0.3491 

radians.  This  problem  can  now  be  formulated  in  the  standard  form  given 
in  Chapter  4;  i.e.,  minimize 


Min  t|»0  =  b9 


(5.3) 


subject  to  discretized  design  constraints 


*i±s  V 


-  b9  <  0  ,  i  =  1, 


(5.4) 


The  constraint  that  the  input  link  b^  is  a  crank  is  imposed  in  the  form 
[26] 


+2  E  lr4  ”  b2l  "  b3  +  bj  <  0 


=  bl  +  b3  -  b2  -  r4  <  ° 


(5.5) 

(5.6) 


Nonnegativity  constraints  on  b  and  b  are  written  in  the  form 

7  8 
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\  =  -b?  <  0 

(5.7) 

%  5  -b8  *  ° 

(5.8) 

The  error  bound  constraint  in  Eq .  5.4  is  imposed  over  a  grid  of 
t  =  19  equally  spaced  points  on  the  range  of  a.  The  driving  kinematic 
constraint  for  this  problem  is  given  as 

<f>l  =  °2  “  ^b5  +  b4  +  al^  =  0  >  i  =  1,  ...  ,19 


5.1.3  Numerical  Results 

5. 1.3.1  Verification  of  Design  Sensitivity  Analysis 

It  is  necessary  to  numerically  verify  that  the  technique  developed 
in  Chapter  4  has  been  correctly  coded.  The  procedure  used  is  briefly 
explained  here. 

The  problem  is  set  up  through  user  supplied  routines  and  input 
data.  The  code  is  allowed  to  run  for  2  design  iterations.  On  the  basis 
of  the  design  sensitivity  vector  calculated  and  the  change  in  design  6b 
obtained  in  the  first  iteration,  changes  in  the  constraint  can  be  pre¬ 
dicted.  This  predicted  value  can  be  compared  with  the  actual  change 
obtained  when  the  constraint  functions  are  evaluated  in  the  second  de¬ 
sign  iteration.  If  the  design  sensitivity  analysis  is  valid,  the  pre¬ 
dicted  and  actual  changes  in  constraints  should  agree  within  a  reason¬ 
able  tolerance.  This  procedure  is  now  used  to  verify  the  design  sensi¬ 
tivity  analysis  for  this  problem.  Consider  grid  point  19,  where  the 
parametric  constraint  of  Eq.  5.4  is  violated  during  a  design  iteration. 
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T 

£1>19  =  [-0.67893,  -.04693,  0.06083,  2.30624,  2.85158,  -1.0, 

0.99853,  0.34923,  -1.0] 

T  r 

The  change  in  design  at  this  iteration  is  6b  =  [0.01663,  0.00115, 
-0.00149,  -0.05650,  -0.06985,  0.02450,  -0.02446,  -0.00853,  -0.2096]. 
The  predicted  change  in  the  constraint  is  thus  given  as 

T 

Aip1  * 19  =  l1*19  6b  =  -0.1832 
P 

The  actual  change  in  the  constraint  function  is  evaluated  as 
1  19 

=  -0.1709 

The  difference  between  the  predicted  and  actual  change  is  6.7%.  The 
design  sensitivity  analysis  is  thus  considered  to  be  valid. 

5. 1.3.2  Optimization  Results 

The  results  obtained  from  the  optimization  procedure  are  as 
presented  in  Table  5.1.  None  of  the  nonparametric  constraints  was 
active  at  the  optimum.  As  given  in  Table  5.1,  there  are  5  critical 
points  on  the  grid  of  the  input  parameter  at  the  optimum  design,  where 
the  error  in  function  generation  is  a  relative  maximum.  The  maximum 
error  obtained  in  Ref.  26  for  this  problem  was  about  0.0019,  which  is 
greater  than  that  obtained  by  the  present  technique. 
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5.2  Example  2:  Kinematic  Synthesis  of  a 
Rigid  Body  Guidance  Mechanism 


5.2.1  Problem  Description 

The  mechanism  shown  in  Fig.  5.2  is  to  guide  the  book-carriage  of  a 
library  book-stacking  vehicle  [28].  It  is  required  that  point  P 
describe  a  path  given  in  Table  5.2.  In  addition  it  is  required  that 
body  8,  the  carrier,  assume  angular  orientation  at  specific  input  grid 
points,  given  in  Table  5.2.  The  mechanism  is  driven  by  the  rotating 
crank,  which  is  body  4. 

Design  variables  b  ,  b  ,  and  b  radially  locate  the  3  joints  of  the 
12  3 

ternary  base  link.  The  angular  orientation  of  these  radial  lines  is 

specified  by  the  three  angles  3  =73.5°,  3  *  61.5°,  and  3  =71.5°. 

1  2  3 

Design  variables  b  to  b  are  the  lengths  of  links  4  to  9,  respectively. 
4  9 

The  length  of  link  10  is  taken  as  10  units.  Other  angles  specified  in 

Fig.  5.2  are  given  as  3  -  “50.0°,  3  -  -63.0°,  and  3  =  -103.0871°. 

4  5  6 

5.2.2  Problem  Formulation 

The  global  x  and  y  coordinates  of  the  coupler  point  P  can  be 
written  in  terms  of  design  and  state  variables  as 


R  =  {x0  +  cos  0o(-75 - A  cos  3y  )  -  A  sin  0O  sin  3,  }  I 

p  o  o  l  4  o  4 


+  (yg  +  Sin  -  a  cos  3^)  +  A  cos  0g  sin  3^ }  J  (5.9) 


Equation  5.9  can  be  symbolically  written  as 


R  =  (R)  I  +  (R  )  J 

p  p  p 

x  y 


(5.10) 
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Figure  5.2  Rigid-Body  Guidance  Problem 
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where  Rp  and  Rp  represent  the  x  and  y  components  of  the  position 
rx  *y 

vector  of  point  P. 


An  artificial  design  variable  b  is  now  introduced  to  impose  the 

10 


following  constraints: 


<Rp  >  "  (Px}d 


cb  ,1-1,  •••,7 

1  10  »  j 


<Ep  >  -  (V<i  1  ‘  C2b10  •  3  '  1 . 7 


<V„  -  4  c,b,„  ,  j  -  1 . 7 


8'g  x  8'd 


3  10 


where  (p  )  and  (p  )  are  the  desired  x  and  y  coordinates  of  point  p, 
x  d  yd 

(0  )  is  the  desired  orientation  of  link  8  and  (0  )  is  the  generated 
8  d  8  g 

orientation  of  link  8. 


In  standard  form,  this  problem  is  formulated  as;  minimize 


*o  5  b10  (5al) 

subject  to  discretized  design  constraints 

(5.12) 

(5.13) 

(5.14) 

The  driving  kinematic  constraint  for  this  problem  is  given  as: 

,d  _  Qd  i  .  . 

(j>l  =  @1  -  <*1  =  0  ,  i=l,  . .  .  ,  7 


5 


4  s 


<RP  )  ■  i  -  cibio  *  0  ■  i  •  1 . 7 


X 


(Rp  ^  ^Py^d  j  c2b10  <  0  *  3  ****  7 


4  - 


<Vg  "  (Vd 


j  c3b10  <  0  »  j  1,  . . . ,  7 


5.2.3  Numerical  Results 


5.2.3. 1  Verification  of  Design  Sensitivity  Analysis 

The  procedure  outlined  in  Section  5. 1.3.1  is  used  here  to  verify 
the  design  sensitivity  analysis.  The  design  sensitivity  vectors  at  a 
certain  iteration  for  the  three  error  constraints  were  obtained  as 
T 

jt1*1  =  [-0.44247,  -0.05346,  0.77940,  0.44687,  0.01196,  -0.06584 

0.16044,  -0.18837,  -1.19525,  -1.0  ] 

T 

£2>4  =  [-1.42060,  0.03326,  0.32397,  -0.13204,  0.48655, 

-0.46958,  0.60349,  0.14956,  -0.65320,  -l.O] 

T 

l3>6  =  [0.24504,  -0.15474,  -0.08915,  0.51890,  -0.82560,  0.20605 

-0.23002,  -0.49248,  0.88695,  -1.0  ] 

The  change  in  design  obtained  at  this  iteration  was 

6bT  =  [-0.00368,  0.02838,  -.001735,  0.002504,  0.01337,  0.02730, 

0.00148,  -.00806,  0.003155,  -0.002811  ] 

The  predicted  changes  in  the  constraints  can  thus  be  computed  as, 

T 

Aib1’1  =  A1* 1  6b  =  -0.00096 
P 

T 

Ai|£’4  =  £2’4  6b  =  -0.00060 
T 

A\|^’6  =  £3,6  6T  =  -0.000014 

The  actual  changes  in  these  constraint  functions  were 

Aip1’1  =  (0.003869  -  0.004834)  =  -0.000965 

A 

a/’4  =  (0.0023898  -  0.002992)  =  -0.000604 
A 

AiJ>3,6  =  (0.00006938  -  0.00008909)  =  -0.0000197 
A 
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Comparison  of  the  actual  and  predicted  changes  in  the  constraints 
shows  that  for  the  first  two  constraints  there  is  a  good  agreement. 
However,  for  the  third  constraint  the  agreement  is  not  as  good.  This 
can  be  attributed  to  numerical  error  in  differencing  very  small  numbers 

5. 2. 3. 2  Optimization  Results 

The  results  from  the  optimization  procedure  are  presented  in  Table 
5.3.  The  first  parametric  constraint  was  active  at  3  critical  points, 
the  second  at  2  critical  points,  and  the  third  at  only  one  critical 
point.  The  maximum  error  obtained  using  the  present  technique  is  20% 
higher  than  that  reported  in  Ref.  28. 

5.3  Example  3  —  Two  degree  of  Freedom 
Function  Generator 

5.3.1  Problem  Description 

The  relationship  u  =  (1  +  v)log  (1  +  w)  is  to  be  mechanized  in 

10 

the  region  0  <  v  <  1,  and  0  <  w  <  1.  This  problem  is  similar  to  the 

numerical  example  given  in  Ref.  6.  The  mechanism  to  be  used  for 

function  generation  is  a  seven-link  mechanism  shown  in  Fig.  5.3.  The 

inputs  v  and  w  are  the  displacements  of  bodies  5  and  6,  respectively. 

These  displacements  are  measured  from  reference  positions,  measured 

along  the  global  axis,  which  are  taken  as  design  variables  b  and  b  • 

1  2 

The  output  u  of  the  mechanism,  is  the  displacement  of  body  2  relative 
to  the  origin  of  the  global  coordinate  system.  This  displacement  is 
measured  along  the  line  of  translation  that  makes  an  angle  y  -  27.881° 
with  respect  to  the  global  x-axis.  the  lengths  of  links  3,  4,  6  and  7 
are  design  variables  b  to  b  .  The  angle  between  links  4  and  6  is 


Table  5.3  Results  for  Example 
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Computing  time  on  IBM  370/168  approximately  3.5  Secs/ Iteration 
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design  variable  b  .  The  other  parameters  related  to  this  problem  are 
3 

as  indicated  in  the  Fig.  5.3. 


5.3.2  Problem  Formulation 

The  displacement  of  the  output  slider,  body  2,  along  the  line  of 
translation  can  be  written  in  terms  of  the  generalized  coordinates  of 
body  2 ,  as 


U  = 


cos  3 


(5.15) 


1 


where  3  =27.881°. 

1 

An  artificial  design  variable  b  is  introduced  such  that 

8 


u  -  u,  .  ,  <  bQ 

g  d  j,k  8 


j  =  1 . 4 

1C  =  1  |  •  •  •  ,  A 


where  (u  )  represents  the  generated  value  of  u,  for  specific  values 
of  the  input  parameters  ocj  and  and  (u^)^  ^  the  desired  value  of 
u,  for  the  same  values  of  the  input  parameters  ct|  and  o^.  The 
function  to  be  generated  in  this  example, 

(ud}j,K  =  (1  +  <4)  log10(1  + 

This  problem  can  be  formulated  in  standard  form  as;  minimize 


*0  "  b8 


(5.16) 


Subject  to  discretized  design  constraints 


.i,1’  k  = 


- (i  +  <4>  ios10o +  “5> 


-  b8  <  °, 


j ,k  =  1,  ...  ,4  (5.17) 
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The  error  bound  constraint  is  imposed  on  a  grid  of  four  equally  spaced 

points,  Fig.  5.4,  on  each  of  the  input  variables  a  and  a  .  Thus  a 

1  2 

total  of  16  grid  points  are  obtained.  No  other  design  constraints  are 
imposed  in  this  problem*  The  driving  kinematic  constraints  for  this 
problem  are  given  as 

<j)^  -  x —  (b^  +  aj)  “0  ,  j  =  1  >  •  •  •  ,4 

d  k 

(f>2  =  xg  -  (b2  +  a2)  =  0  ,  k  =  1 , ...  ,4 

5.3.3  Numerical  Results 

5.3.3. 1  Verification  of  Design  Sensitivity  Analysis 

The  design  sensitivity  of  the  upper  bound  constraint,  Eq.  5.17,  at 
iteration  no.  11  and  grid  point  j  =  1,  k  =  4  was  obtained  as 
T 

A1*1*4  =  [0.15546,  0.23616,  -0.57531,  -1.45199,  1.29105,  0.06958, 

-0.32795,  -l.O] 

The  change  6b  in  design  was 

6bT  =  [-0.01027,  0.00197,  0.01753,  0.00938,  0.01017,  0.01178, 
-0.00640,  -0.00070] 

The  predicted  change  in  constraint  is  thus  computed  as 

T 

A*1*1*4  =  A1*1*4  6b  =  -0.00805 
P 

The  actual  change  in  constraint  was  obtained  as 


A^A* 1 *4  =  °*07259  -  0.08055  =  -0.00796 
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The  difference  between  the  actual  and  predicted  change  in  the  con¬ 
straint  is  1.0%,  which  verifies  that  the  design  sensitivity  analysis 
for  upper  bound  constraint  of  Eq.  5.17  is  valid. 

5.3.3. 2  Optimization  Results 

The  results  of  the  optimization  procedure  for  this  problem  are  as 
presented  in  Table  5.4.  The  maximum  error  in  function  generation 
obtained  from  the  present  procedure  is  about  25%  higher  than  that 
obtained  in  Ref.  6.  This  can  be  attributed  to  2  differences  in  the 
formulation.  First,  the  design  variables  used  in  the  present  problem 
and  Ref.  6  are  not  the  same.  Second,  the  grid  points  in  the  present 
formulation  and  Ref.  6  are  not  identical.  Location  of  grid  points  is 
known  to  have  a  very  significant  effect  on  the  error  in  function 
generation  [14]. 


5.4  Example  4  -  Stress  Constrained  Design 
of  a  Four-Bar  Mechanism 

5.4.1  Problem  Description 

The  four-bar  mechanism  shown  in  Fig.  5.5  has  its  input  crank 
rotating  at  a  constant  angular  velocity  of  300  rpm.  This  mechanism  is 
to  be  designed  for  minimum  mass,  requiring  that  bending  stresses  in 
mobile  links  2,  3,  and  4  do  not  exceed  an  allowable  stress.  The  design 
variables  in  the  mechanism  are  the  circular  cross-sectional  areas  of 
the  three  mobile  links.  The  link  lengths  are  kept  fixed  at  the  values 
specified  in  Table  5.5.  The  mass  density  and  allowable  stress  in  the 
links  are  as  given  in  Table  5.5.  Specifications  for  this  problem  are 
the  same  as  those  of  example  1  in  Ref .  30  and  example  1  in  Ref.  29. 


Figure  5.5  Minimum  Weight  Design  of  Four-Bar  Mechanism 
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5.4.2  Problem  Formulation 

As  shown  in  Chapter  3,  the  Lagrange  multipliers  obtained  as  a 

solution  of  the  force  equation  are  directly  related  to  reaction  forces 

in  the  joints.  As  shown  in  Fig.  5.6,  reactions  at  the  ends  of  a  link 

could  thus  be  represented  by  the  appropriate  Lagrange  multipliers. 

However,  to  compute  bending  moments,  the  end  reactions  are  required  to 

be  expressed  in  a  local  coordinate  system.  As  shown  in  Fig.  5.6,  the 

p  ,  U  system  of  forces  at  a  joint  must  be  transformed  to  the  P  , 
k  k+1  k 

P  system. 
k+1 

The  external  load  on  the  beam  depends  upon  the  normal  acceler¬ 
ation.  To  compute  the  distribution  of  normal  acceleration  along  the 
length  of  the  link,  a  cross  section  of  the  link  at  distance  v  along  the 
x  axis  is  considered.  The  position  vector  for  point  c  can  be  written 
as 


R  =  x  I  +  y  J  +  v(cos  0.  I  +  sin  0  J)  (5.18) 

Ci  111 

where  x  ,  y^  and  0^  are  the  three  generalized  coordinates  locating  the 
body  in  the  plane  and  I  and  7  are  unit  vectors  in  the  global  X  and  Y 
directions. 

The  absolute  acceleration  of  point  c  can  thus  be  written  as 
R  =  (x.  -  v  cos  6.  0^  -  v  sin  0.  0.)  I  + 

c±  l  l  i  ii 

•*  *2  ••  — 

(yi  -  v  sin  0±  0.^  +  v  cos  0^  )  J 

Transforming  the  global  unit  vectors  to  the  local  system, 
acceleration  normal  to  the  length  of  the  link  can  be  written  as 


H-X+l 
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Figure  5.6  Loading  System  in  Problem 
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•  •  *• 

a  =  (-  x.  sin  0.  +  y.  cos  0.)  +  v  0.  (5 • 19) 

y.  l  l  l  i  x  ' 

where  x^,  y^  and  0^  are  the  generalized  accelerations  of  body  i.  From 
Eq.  5.19,  it  can  be  seen  that  the  normal  acceleration  varies  linearly 


along  the  length  of  the  link. 


The  distributed  loading  on  the  beam  can  thus  be  written  as 


f.  =  {(-x.  sin  0.  +  y.  cos  0.)  +  v  0.  }  (5.20) 

i 


where  ra  is  mass  per  unit  length  of  link  i. 
I 

l 

Equation  5.20  can  be  written  as 


q  =  -m*.  a*.  +  v  ei 

i  i 

where  a^  =  sin  0^  +  y^  cos  0^. 


(5.21) 


Since  elementary  beam  theory  is  used  to  compute  bending  stresses, 

it  is  necessary  to  first  determine  the  point  of  maximum  bending  moment. 

This  is  the  point  at  which  the  shear  force  changes  sign.  The  shear 

force  at  any  point  D,  at  a  distance  X  from  0,  can  be  written  as 

0 


F 

s 


*  P 


k+1 


+ 


+  v  0. )  d  v 
i 


=  P 


k+1 


t) 


2 


(5.22) 


At  the  point  of  zero  shear  force,  the  right-hand  side  of  Eq.  5.22 
will  be  zero,  giving  rise  to  the  condition 
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a*A 


0 

-4  *1 +  a*.  xo +  H 


2  " 
Z.  0. 
i  l 

8 


k+1 


)  =  0 


(5.23) 


Solving  for  X  from  Eq.  5.23, 
0 


X0  = 


20. 

1 


a*A 

l 


2  ’• 
\  ei 
8 


k+1 
m  n 


(5.24) 


The  maximum  bending  moment  can  thus  be  written  as 


Z. 


A  pk+l  (2  +  xo)+mz  fg/  (aZ.  +  V  6i  Kn  -  x0)  dv 

*■  &/  2  i 


(5.25) 


Integration  of  Eq.  5.25  gives 


Z  T  X?  Z*  Z.X„ 

\  =pk+J-2  +  V  +  “ik[-T  (-ji +  -+£-)] 


X  z3  X  z2  1 
+  0  r — 2.  +  -JL  +  1 

i  L  6  24  +  8  J I 


(5.26) 


where  X^  is  given  by  Eq.  5.24. 

Equation  5.26  is  valid  only  for  links  3  and  4.  The  crank,  link  2, 

is  rotating  at  constant  angular  velocity  and  has  no  normal 

acceleration.  The  maximum  moment  on  this  link  is  the  torque  required 

* 

to  drive  the  mechanism.  It  occurs  at  i)  =  -  and  is  given  as 


M  =  u 
02  12 


(5.27) 


Since  the  links  are  of  circular  cross  section,  the  area  moments  of 
inertia  can  be  written  in  terms  of  the  design  variables  as 
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(5.28) 


where  2  <  i  <  4  is  the  link  number.  The  distance  of  the  extreme  fiber 
from  the  neutral  axis  is  given  as 


The  absolute  values  of  the  maximum  bending  stresses,  from 
elementary  beam  theory  [31],  can  now  be  written  as 


(Vi  = 


,3/2 

Vl 


2  <  i  <  4 


(5.29) 


where  Mq  is  given  by  Eq .  5.27  for  i  =  2  and  by  Eq.  5.26  for  i  =  3  and  4 

The  optimization  problem  can  be  written  in  standard  form  as; 
minimize 


*0  “  P2  %2  bl  +  p3  *3  b2  +  P4  \  b3 


(5.30) 


where  is  mass  density  of  link  i,  subject  to  discretized  design 
constraints 


(<4)k+l 


a 

a 


1  <  0  ,  j  = 


1, 


T 


1  <  K  <  3 


(5.31) 


where  is  allowable  stress  and  x  is  the  total  number  of  grid  points 
For  this  problem  19  grid  points  were  considered.  The  driving 
constraint  for  this  problem  is  given  as 

*1  =  02  "  al  =  0  »  j  =  1,...,19 


5.4.3  Numerical  Results 


5.4.3. 1  Verification  of  Design  Sensitivity  Analysis 

The  design  sensitivities  of  the  stress  constraints  at  a  certain 
iteration  were  obtained  as 

T 

A1’17  =  [-1801.18,  1928.83,  853.239]1 

T 

A2’15  =  [  0.0  ,-3481.67,  0.0  ]T 

T 

£3’15  -  [  0.0  ,  0.0  ,+2051.6  ]T 

The  change  in  design  at  this  iteration  was 

6bT  =  [2.022xl0”5,  1.76xl0"5,  -8.626xlo"6] 

The  predicted  changes  in  the  constraints  can  thus  be  computed  as 

T 

A^’17  =  A1’17  6b  =  -0.00983 
T 

A^>15  =  l2>15  6b  =  -0.0613 
T 

A^’15  =  A3’15  6b  =  -0.01739 

The  actual  changes  in  the  constraint  function  were  obtined  as 

a4’17  =  -o*°097 

A1^>15  =  -0.06002 
Aip3’15  =  -0.0169 

Comparison  of  the  actual  and  predicted  changes  shows  that  the  two 
sets  of  data  agree  to  within  2.5%.  The  design  sensitivity  analysis  can 


thus  be  considered  valid 
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5. 4.3. 2  Optimization  Results 

Results  of  the  optimization  procedure  are  presented  in  Table  5.6. 
Results  obtained  in  Refs.  29  and  30  for  this  problem  are  also  given. 

As  can  be  seen,  there  is  a  substantial  difference  between  the  designs 
obtained.  This  difference  is  attributed  to  the  fact  that  Refs.  29  and 
30  use  different  methods  for  stress  analysis  than  used  in  the  present 
formulation.  The  constraint  functions  are,  therefore,  not  identical. 

5.5  Example  5  -  Dynamic  Balancing  of 
a  Four-Bar  Mechanism 

5.5.1  Problem  Description 

It  is  desirable  to  make  the  shaking  forces  of  a  mechanism  vanish 
by  balancing  the  mechanism.  However,  doing  so  introduces  additional 
difficulties.  For  example,  the  RMS  bearing  forces  could  increase  by 
as  much  as  100  percent  [4].  In  the  present  example,  the  four-bar 
mechanism  shown  in  Fig.  5.7  is  considered.  A  trade-off  is  resorted  to, 
in  which  the  RMS  shaking  forces  are  minimized,  while  limiting  increase 
in  the  ground -bearing  forces.  Minimization  of  RMS  shaking  force  is 
accomplished  by  modifying  the  inertial  properties  of  the  output,  link 
by  adding  a  circular  counterweight  to  it,  as  shown  in  Fig.  5.8.  The 
mass  and  location  of  this  counterweight  are  design  variables. 

Data  for  the  mechanism  considered  in  this  problem  are  given  in 
Table  5.7.  This  problem  is  the  same  as  the  single  counterweight 
optimization  example  given  in  Ref.  4. 


Table  5.6  Results  for  Example  4 


Design  Variables 

Mass 

(kg) 

bi(m2) 

b2(m2) 

b3(m2) 

Starting 

Design 

1x10-1 

1x10-1 

1x10-1 

546.26 

Optimal  Design 
30th  iter 

1.03x10-3 

0.621x10-3 

0.1169x10-3 

2.686 

Results  from 
Reference  29 

0.401x10-3 

0.385x10-3 

0.385x10-3 

2.134 

Results  from 
Reference  30 

0.953x10-3 

0.488x10-3 

0.309x10-3 

2.703 

At  the  final  design  |  | 6b  |  |  =  2.5  x  10 


Approximate  CPU  time  (IBM  370/168)  =  1.25  Sec/Iteration 
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Table  5.7  Data  for  Example  5 


LINK 

12  3  4 

Length  (m) 

7.62x10-2 

2.54x10-2 

1.016x10-1 

7.62x10-2 

Mass  (kg) 

— 

4.588x10-2 

1.084x10-1 

6.602x10-2 

Moment  of 
inertia  @ 
C.G  kg-m2 


3.208x10-4 


6.779x10-5 
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5.5.2  Problem  Formulation 

The  RMS  shaking  force  for  the  four-bar  mechanism  being  considered 
here  is  given  as  [4] 

(FM/g)rms  =  4)  ^F14X  +  F34X^  +  ^F14Y  +  F34Y^  Jd02 

(5.32) 

where  F  and  F  are  the  ground-bearing  forces  at  joint  1,  in  the  X 
14X  14Y 

and  Y  directions,  respectively;  F  and  F  are  the  ground- bearing 

34X  34Y 

forces  at  joint  4,  in  the  X  and  Y  directions,  respectively;  and  0  is 

2 

the  angular  orientation  of  the  input  crank. 

The  RMS  ground-bearing  forces  are  given  as  [4], 


(5.33) 


(5.34) 


where  (f  )  and  (f  )  are  the  RMS  ground -bearing  forces  at  joints 
1  RMS  4  RMS 

1  and  4,  respectively. 

There  is  a  one  to  one  correspondence  between  the  design  variables 
used  here  and  those  used  in  Ref.  4.  Location  of  the  center  of  mass  of 
the  combined  link,  in  terms  of  design  variables  and  given  data,  can  be 
written,  using  relationships  given  in  Ref.  4,  as 

°4  ■  p4  ‘  ToTTT  (m4p4  -  blb2>  (5‘35) 

K  +  bi  J 


and 
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(m4  +  bi ) 


(5.36) 


where  c  and  d  define  the  location  of  the  combined  center  of  mass,  b 
4  4  1 

is  the  mass  of  the  counterweight,  b  and  b  are  the  distances  shown  in 

0  2  3  o  _ 

Fig.  5.8,  m  is  mass  of  the  original  output  link,  and  p  is  given  data, 
4  4 

as  defined  in  Fig.  5.8. 

Acceleration  of  point  A  relative  to  point  0  can  be  computed  by  the 
relation 


aA0  2  (rA()) 
dt 


(5.37) 


where  r  is  the  position  vector  of  point  A,  relative  to  point  0,  given 
AO 


rA0  =  (c4  cos04  -  d4sin04)I  +  (c4sin04  +  d4cos04)j  (5.38) 


Substituting  r  from  Eq .  5.38  into  Eq .  5.37  gives 
AO 


aAo  =  Hca  cos 


e/.  "  d/.sin6/.  K®/J2  "  (c4  sin04  "  d4cos64)04}1 


{-(c4  sin04  +  d4cos04  )(04  )2  +  (c4  cos04  -  d4sin04)0‘4KJ 


.  w  o  ,  «/  U4.UW,  7W/' 

4  4  4  4 J  4 


(5.39) 


The  absolute  acceleration  of  point  A  can  thus  be  written  as 


aA  a0  +  aA0 


(5.40) 


D'Alembert's  force  on  the  output  link  can  now  be  written  as 


F4  ~  ”(m4  +  bi)  aA 


(5.41) 


Computation  of  D'Alembert  moments  requires  the  combined  moment  of 
inertia  of  the  original  link  and  the  counterweight.  Since  the 
counterweight  is  circular,  it's  moment  of  inertia  about  point  B  is 
given  as 

I*  =  b^b2  +  b2)/2  (5.42) 

Using  the  parallel  axis  theorem  [23],  the  combined  moment  of  inertia  of 
the  original  link  and  the  counterweight  can  be  written  as 

IA-I0+“4(l^l)2  +  I‘  +  bl<l;ABl)2  «•«> 

where  I  is  the  combined  moment  of  inertia  about  point  A  and  I  is  the 
A  0 

moment  of  inertia  of  the  original  link  about  point  0.  Since  moment  of 

inertia  of  the  original  link  given  in  Ref.  4  is  about  point  C,  the 

parallel  axis  theorem  must  be  used  to  write  I  in  Eq.  5.43  in  terms 

0 

of  the  given  data, 

IC  =  I0  +  m4  (p4^2  (5.44) 

Substituting  for  I  in  Eq .  5.42  from  Eq .  5.44  gives 
0 

\  ■  xc  “  “4^4)  +  m4(c4  +  d4> 

+  I*  +  b1  [(p°  -  c4  +  b2)2  +  (b3  -  d4)2]  (5.45) 


The  D’Alemberts  moment  can  now  be  written  as 


100 

Forces  in  the  ground  bearings,  used  inEqs.  5.32  to  5.34,  can  be 
related,  from  the  results  of  Section  3.3.1,  to  the  Lagrange  multipliers 
by  the  following  relationships 


14X  ~ 

"P4 

14Y  “ 

"y5 

34X  ” 

yio 

34Y  ~ 

P11 

(5.47) 


Equations  5.32  to  5.34  can  thus  be  rewritten  as 


-  J-k  ?Vh*  +  V*  +  <-^5  +  |1n)2]  d62  <5-48) 


I y  +  ^  ]  d9„ 

4  5  J  2 


I  i  2ir  : 

l^RMS  =  “2i  ^  ^ 

f  211 

4)RMS  =  {  ^0  +  ^  d02 


(5.49) 


(5.50) 


The  optimization  problem  can  now  be  stated  in  the  standard  form 
of  minimizing 


*0  "  ^m/gVs 


(5.51) 


subject  to  constraints 


tl  $  (VbMS  .  !  <  o 


(5.52) 


(F  ) 

,2  _  4  RMS 

*  =  —5 - 


1  <  0 


(5.53) 
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where  (Fm/q)^,  (F^g,  and  (F4)RMS  are  given  by  Eqs.  5.48,  5.49, 

and  5.50,  respectively,  and  F  and  F  are  the  upper  bounds  on  the 

1  4 

bearing  forces  at  revolute  joints  1  and  4,  respectively. 

The  driving  kinematic  constraint  for  this  problem  is 

<j)j  =  0  2  —  cc|  =  0  ,  j  =  1, ...  ,15 


The  rationale  used  to  compute  the  values  of  F  and  F  is  explained 

1  4 

in  detail  in  Ref.  4.  Briefly,  to  arrive  at  these  numbers  the  RMS 

ground  bearing  forces  for  the  unbalanced  mechanism  and  the  fully  force 

balanced  mechanism  are  computed. 

Since  the  RMS  bearing  forces  for  the  fully  force  balanced 

mechanism  are  higher  than  those  for  the  unbalanced  mechanism,  the 

values  of  F  and  F  are  chosen  to  lie  midway  between  the  two  extreme 
1  4 

values.  For  the  present,  these  values  are  taken  as  F^  =  9.136  *  10~3  N 

and  F  =  6.427  x  10"3  N. 

4 

The  cost  function  and  design  constraints  considered  in  example 
problems  thus  far  do  not  involve  integrated  quantities  such  as  those 
appearing  in  Eqs.  5.48  to  5.50.  Computation  of  design  sensitivities  of 


these  functions  can  also  be  handled  by  the  formulation  developed  here. 
Design  sensitivities  of  the  integral  cost  and  constraint  functions  can 
similarly  be  written  by  taking  the  derivatives  with  respect  to  design 


of  Eqs.  5.51  to  5.53  [32],  yielding 
T  dak 


l°‘  =  ° 


db  2it(Fm/G^EMS  0 


2  or  d  y.  d  p.  n 

I  +  “db +  “dr^ 


d  dp.. 

+  (-u5  +  Uu>(--db  +-db")]  d02 


(5.54) 


102 


£1  _  dip  _ 

4  "  dF“  “ 


If,  2  IT 

F 1  |  2tt(F  )  ^  t  H 

L  v  i'rms  o 


2  IT  d  M  dp 

‘4  ~~db  +  y5  ~db d62 


(5.55) 


/.i i_r  r_i _ **,  %,  +  B  ibild9]  ..... 

db  F4  )_2tt(F4)rms  '  Ly10  db  +  yll  db  J  d02j  (5,56) 


Numerical  integration  of  Eqs.  5.54  to  5.56  requires  evaluation 

of  the  integrand  at  specific  points  in  the  interval  of  integration. 

Computation  of  the  derivatives  of  Lagrange’s  multipliers  appearing  in 

the  integrands  of  Eqs.  5.54  to  5.56  can  be  routinely  done  in  the 

present  formulation.  The  computer  code  is  forced  to  treat  y  ,  y  ,  y  , 

4  5  10 

and  y  as  constraint  functions,  purely  for  the  purpose  of  computing 

the  design  sensitivity.  The  design  sensitivity  of  these  functions  is 

evaluated  at  points  on  the  0^  grid  corresponding  to  the  nodes  of  the 

15  point  Gauss-Legendre  integration  formula  [33].  Since  constraints 
1  2 

^  and  ^  do  not  depend  on  the  input  parameter,  they  are  nonparametric 
constraints. 


5.5.3  Numerical  Results 

5.5.3. 1  Verification  of  Design  Sensitivity  Analysis 

The  design  sensitivity  vectors  of  the  cost  and  constraint 
functions,  at  a  certain  design  iteration,  were  obtained  as 
T 

1°  =  [-.000799,  .02880,  .004828] 

1T 

r  =  [4.1542,  -12.7106,  3.1470] 

9T 

r  =  [3.9492,  -14.846,  -.3122] 
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The  change  in  design  in  this  iteration  was 

6b  =  [-0.000885,  0.001334,  0.000373] 

The  predicted  change  in  cost  is  thus  computed  as 
T 

Aip°  =9°  6b  =  0.0000409 
P 

Predicted  changes  in  constraints  are  computed  as 

Aip1  =  -.01946 
P 

A*2  =  -.0234 
P 

The  actual  changes  in  the  cost  and  constraints  were 
Ai|>°  =  0.000057 
=  -.011355 
A<|>2  =  -.01406 

Comparing  the  predicted  and  actual  changes  in  cost  and  constraint 
functions,  it  can  be  seen  that  the  two  sets  of  data  no  not  agree  well. 
This  can  be  attributed  to  the  highly  nonlinear  nature  of  the  cost  and 
constraint  functions.  This  can  be  seen  by  computing  the  predicted 
change  in  constraint,  using  an  average  design  sensitivity  vector.  For 
example,  consider  constraint  1,  for  which  the  design  sensitivity  vector 
averaged  over  the  two  design  iterations  is 

iT 

[V  -  [4.1145,  -12.4848,  3.1727  ] 

L  JAVE  L 

The  predicted  change  in  the  constraint,  on  the  basis  of  the  averaged 
design  sensitivity  vector,  is 


104 


KJave  ■  --01183 

As  can  be  seen,  this  figure  agrees  closely  with  the  actual  change  in 
the  constraint. 

5. 5. 3. 2  Optimization  Results 

Results  obtained  from  the  optimization  procedure  are  presented 
in  Table  5.8.  For  the  purpose  of  comparison,  the  design  variables  used 
in  Ref.  4  are  converted  to  those  used  in  the  present  formulation.  The 
ratio  of  the  RMS  shaking  force  for  the  partially  balanced  mechanism  to 
that  for  the  unbalanced  mechanism  is  given  by  r  .  As  can  be  seen  from 
Table  5.8,  a  substantial  reduction  in  RMS  shaking  force  is  obtained  by 
using  the  present  formulation.  Results  obtained  from  the  present 
formulation  and  those  in  Ref.  4  differ  substantially.  This  difference 
is  attributed  to  the  fact  that  the  present  formulation  imposes  the  RMS 
ground  bearing  force  constraints  as  inequality  constraints  and  Ref.  4 
considers  these  to  be  equality  constraints.  Significant  improvement  in 
design  from  the  present  formulation  over  that  from  Ref.  4  is  thus  to  be 
expected  and  is  evident  from  Table  5.8. 


05 


Table  5.8  Results  in  Example  5 


Design  Variables  in 
present  paper 

b^kg) 

b2(m) 

b^(m) 

Starting  Design 

0.30 

0.0 

0.0 

— 

Optimal  Design, 
15th  Iter 

0.203 

0.040 

0.0098 

0.25 

Optimal  Design 
Ref.  4 

0.207 

0.0183 

-0.017 

0.69 

At  the  final  design  j  | 6b  | |  =  1.6  x  10 
Approximate  computing  time  IBM  370/168  =1.0  Secs/iteration 
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CHAPTER  6 

CONCLUSIONS  AND  RECOMMENDATIONS 
FOR  FURTHER  RESEARCH 

6*1  Conclusions 

The  five  examples  presented  in  Chapter  5  were  studied  to  evaluate 
versatility  and  effectiveness  of  the  technique  developed  in  the 
previous  chapters.  In  light  of  the  results  obtained  with  these 
examples,  the  following  conclusions  may  be  drawn: 

1.  The  constrained  multi-element  technique  developed  in 
Chapters  2  and  3  for  modeling  kinematic  systems  is  general 
enough  to  represent  a  wide  variety  of  planar  mechanisms, 

2.  The  adjoint  variable  technique  has  been  successfully  applied 
for  design  sensitivity  analaysis  of  planar  mechanisms, 
considering  both  kinematic  and  force  response  of  mechanisms. 

3.  Due  to  the  generality  of  the  constrained  multi-element 
technique  for  modeling  and  the  adjoint  variable  technique  for 
design  sensitivity  analysis,  the  synthesis  and  design 
technique  developed  in  this  research  is  very  general  and 
versatile.  This  is  evident  from  the  variety  of  examples 
given  in  Chapter  5. 

6*2  Recommendations  for  Further  Research 

In  the  preceding  section  it  was  pointed  out  that  the  state  space 
technique  developed  in  this  research  has  demonstrated  its  potential 


as  a  versatile  and  powerful  tool  for  kinematic  synthesis  and  design. 
However,  to  fully  realize  this  potential,  it  is  recommended  that  the 
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following  areas  be  further  investigated: 

1.  The  gradient  projection  algorithm  stated  in  Chapter  5  was 
used  as  the  optimization  algorithm  in  this  research.  This 
algorithm,  in  its  present  form,  is  insensitive  to 
nonlinearity  of  the  state  equations  and  in  the  cost  and 
constraint  functions  encountered  in  kinematic  synthesis  and 
design  problems.  In  addition  to  its  insensitivity  to 
nonlinearities,  this  optimization  algorithm  often  produces 
designs  for  which  the  mechanism  locks  up  in  some  position. 

This  is  primarily  due  to  the  algorithm's  not  having  any 
information  about  the  kinematics  of  the  mechanism. 

It  is  thus  desirable  to  investigate  alternate  optimization 
algorithms  that  have  the  following  desirable  characteristics: 
(i)  The  algorithm  should  be  sensitive  to  nonlinearities 
arising  in  mechanism  synthesis  and  design  problems. 

(ii)  The  algorithm  should  produce  design  changes  that  are 

compatible  with  kinematic  requirements  of  the  mechanism. 

2.  The  gradient  projection  algorithm,  like  many  other  opti¬ 
mization  algorithms,  does  not  guarantee  a  globally  optimal 
design.  To  be  assured  that  the  optimal  design  produced  by  the 
optimization  technique  is  globally  optimal,  it  is  necessary  to 
start  the  procedure  with  different  initial  designs.  In  most 
mechanism  synthesis  problems  a  substantial  effort  has  to  be 
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put  into  producing  acceptable  initial  designs.  It  is  thus 
desirable  to  investigate  techniques  that  guarantee  conver¬ 
gence  to  globally  optimal  designs, 

3.  The  computer  code  used  to  implement  the  theory  developed 
in  the  previous  chapters  was  of  an  experimental  nature.  Some 
of  the  areas  in  which  this  code  could  be  refined  are  as 
follows: 

(i)  Significant  improvements  should  be  made  in  the  solution 
procedure  for  state  equations. 

(ii)  Computation  of  derivatives  of  user-supplied  driving 

kinematic  equations  and  cost  and  constraint  functions 
should  be  carried  out  with  a  symbolic  manipulation 
language . 

(iii)  Implementation  of  the  code  for  execution  in  a 
conversational  mode  would  be  highly  desirable. 

4,  Synthesis  and  design  of  three-dimensional  mechanisms  is 
currently  at  an  elementary  level.  The  theory  developed  here 
for  synthesis  and  design  of  planar  mechanisms  is  extendable  to 
three-dimensional  mechanisms  and  should  be  pursued. 
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APPENDIX  A 

PROOF  OF  NONSINGULARITY  OF  MATRIX  A  IN  EQUATION  4.27 


Matrix  A  defined  in  Eq.  4.27  is  the  coefficient  matrix  on  the 
left-hand  side  of  Eq.  4.26.  A  proof  of  the  existance  of  an  inverse  for 
this  matrix  will  now  be  developed.  The  Jacobian  matrix  (3$/3z)  that 

j 

occupies  a  band  along  the  diagonal  of  matrix  A  is  assumed  to  be 
nonsingular.  This  assumption  is  initially  made  in  Section  2.2.3  and  is 
still  valid  here. 

Consider  Eq .  4.22,  with  perturbation  in  design  6b  =  0, 


3$ 

3z 

Since  the  coefficient  matrix  of  Eq.  A.l  is  nonsingular,  by  theorem  7.2 
of  Ref.  21,  Eq.  A.l  implies  that 

Sz"*  =  0  (A. 2) 

Consider  Eq.  4.23  with  6b  =  0, 


6z 


=  0 


(A.l) 


lz 


r 

“ 

. 

. 

7 

— 

3$  z 

+  9 

3$  a 

J_9z 

3z 

9z 

3a 

I 

6z 


(A. 3) 


Substituting  for  6z^  from  Eq.  A. 2,  Eq.  A. 3  becomes 


3$ 

lz 
L.  Jj 


5z  =0 


(A. 4) 
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By  the  same  argument,  Eq.  A. 4  implies 
6zj  =  0 

Consider  Eq.  4.24  with  6b  =  0  and  Eqs.  A. 2  and  A. 5, 


(A. 5) 


[«! 


6z  =0 


(A. 6) 


As  deduced  previously,  Eq.  A. 6  implies 

6z  =  0 

Consider  Eq.  4.25  with  6b  **  0  and  Eqs.  A. 2,  A. 5  and  A. 7, 

r  n  T  J 

3$  <$y  a  0 

3z 


(A. 7) 


(A. 8) 


The  coefficient  matrix  of  Eq.  A. 8  is  the  transpose  of  the  Jacobian 
matrix  of  Eq.  A.l.  Since  the  Jacobian  is  nonsingular,  its  transpose 
also  has  the  same  property.  Eq.  A. 8  thus  implies 


6  y  =  0 


(A. 9) 


From  the  definition  of  the  composite  state  variable  vector  U  given 

in  Section  4.3, 

_  - 
•  T  rji  ^rji  •••£  iji 

6U  =  6z  6z  6z  6y  (A. 10) 

~  j 

Substituting  from  Eqs.  A. 2,  A. 5,  A. 7  and  A. 9,  Eq.  A. 10  can  be  written 


6  U  -  0 

Now  consider  Eq .  4.27  with  6b  =  0, 


(A. 11) 


j  3  - 

A  6  U  =  0 


(A. 12) 


As  shown  in  the  above,  Eq .  A. 12  implies  that  Eq .  A. 11  holds,  i.e*  , 

j 

6  U  =0 

By  theorem,  7 *2  of  Ref*  21,  this  is  equivalent  to  the  statement  that 
(A^)  is  nonsingular. 
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ac  .  an. .  as. .  an. . 

The  partial  derivatives  _ *j_,  _ ^L,  jjj  and  Jji,  will  take  on 

"5b 

k  k  k  k 

nonzero  values  only  if  design  variable  b  is  related  to  the  lengths  of 

k 

body  i  or  j . 

Similarly  for  translational  joint,  the  following  notation  is 


introduced: 


T1  =  K  COS  0  -  n  sin  0 

ij  i  ij  i 


T2  =  5  sin  0  +  n  cos  0 

ij  i  ij  i 


T3  =  5  cos  0  -  n  sin  0 

ji  j  ji  j 

\  (B.3) 

T4  =  K  sin  0  +  n  cos  0  ( 

ji  3  Ji  j 


T5  =  (T1)(T3)  +  (T2)(T4) 


T6  =  -  (T2)(T3)  +  (T1)(T4) 
T7  =  x  -  x 

i  j 


T8  =  y  -  y 

i  j 
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TB4  =  914  =  J1  sin  6.  +  31  cos  0. 

"55  9b  j  9b  j 

k  k  k 


TB5  =  (TB1) (T3)  +  (T1)(TB3)  +  (TB2)(T4)  +  (T2)(TB4) 


TB6  =  -  (TB2) (T3)  -  (T2)(TB3)  +  (TB1)(T4)  +  (T1)(TB4) 


35..  3^  .  3C  3H 

The  partial  derivatives  ^  ^  ,  j  and  j  in  Eq.  B.4 

"55  "55  "55  155 

k  k  k  k 

will  take  on  non-zero  values  only  if  the  design  variable  b  is  related 

k 

to  location  of  the  translational  joint  connecting  bodies  i  and  j. 


B.2  Derivatives  With  Respect  To  Design  Variables  Only 


Consider  now  each  of  the  two  types  of  joints  separately. 


B.2.1  Revolute  Joint 


From  Eqs.  2.10  the  derivative  of  constraint  equations  with  respect 


to  design  variables  could  be  written  as 
9  4> 

_ *  =  RBI  -  RB3 

9b 

k 
9  <f> 

_ L  =  RB2  -  RB4 

9b 

k 


where  RBI,  RB2,  RB3  and  RB4  are  defined  by  Eq.  B.2 


B.2.2  Translational  Joint 


From  Eqs.  2.13  and  2.16,  the  derivatives  of  the  constraint 


equations  for  the  translational  joint  with  respect  to  design  can  be 


written  as 
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9  (j) 

n  =  (TB1) (U  -U  )  +  (U  -x  )(TB1-TB3)  + 
8b  i  j  i  i 

k 

(TB2) (V.-V.)  +  (V . -y  ) (TB2-TB4) 

l  j  li 


9  (() 

_ ®  =  (TB1 ) (V  -y  )  +  (U  -x  ) (TB4)  - 

9b  3  3  i  i 

k 

(TB2)(Uj-x^)  +  (Vi-yi)(TB3) 


(B.6) 


where  TB1 ,  TB2,  TB3  and  TB4  are  given  by  Eqs.  B.4  and  U  ,  U  ,  V  and 

i  j  i 

V  are  given  by  Eqs.  2.14. 

3 


B.3  Derivatives  With  Respect  To  State  Variables 

B.3.1  First-Order  Partial  Derivatives 

B.3. 1.1  Revolute  Joint 

From  Eq .  2.10,  the  nonzero  first-order  partial  derivatives  of 
equations  of  constraint  for  the  revolute  joint  can  be  written  as 


9  (j> 
x 


9x 

i 


=  1 


9<J> 

x 

90~ 

i 


-  R2 


-1 


9  (j> 


x 


90 

j 


=  R4 


(B.  7) 
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3  (j> 

_z  =  i 
9y 

i 

3  <J> 

_y  -ri 

99 

i 


3  <f> 


Z  =  -  l 


9y. 


9  4> 


Z  =  -  R3 


99 


B. 3«1.2  Translational  Joint 

From  Eqs.  2.13  and  2.16,  the  nonzero  first-order  partial  derivatives 
of  the  equations  of  constraint  for  a  translational  joint  can  be  written 


as 

9  4> 

_ “  =  T1 

3x 

i 


3  $ 

_n  =  T2 

9y 

i 


9<|> 

=  -  (T2)(T7)  +  (T1)(T8)  -  T6 
99 


9* 

n  =  -  T1  >  (B.8) 

9x  ( 

3 


9<j) 

n 

9y~ 

j 


-  T2 
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9<f> 

__H  =  T6 
99 

j 

_ 1  =  -  T5 

99 

i 


80 

j 


B.3.2  Second-Order  Partial  Derivatives 

B.3.2.1  Revolute  Joint 

Differentiating  the  derivatives  in  Eq.  B.7  once  with  respect  to 
the  state  variables  gives  the  following  nonzero  second  derivatives: 


924> 

_ *  =  -  R1 

2 

99 

i 

92  cf> 

_ 5.  =  R3 

2 

99 

j 

924. 

_ y  =  -  R2 

2 

99 

i 

92  <J> 

_ y  =  R4 

2 

99 

j 


(B.9) 
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B«3.2»2  Translational  Joint 

Differentiating  the  derivatives  in  Eq.  B.8  with  respect  to  state 
variables  gives  the  following  nonzero  second  derivatives: 


2 

d  <j> 


n 


9x  90 
i  i 


=  -  T2 


2 

9  (j) 


n 


9y  99 
i  i 


T2 


2 

9  (j> 


n 


99  9x 
i  j 


2 

9  (j> 
n 

99  9y” 

i  j 


2 

9  <(> 
n 

99  99 

i  J 


2 

9  <j> 
n 


90 


=  T2 


=  -  T2 


=  T5 


(B.10) 


=  -  (T1)(T7)  -  (T2)(T8)  +  T5 


2 

9  <j> 

_ R  =  T5 

2 

90 

j 
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B.3.3  Third-Order  Partial  Derivatives 


B.3.3.1  Revolute  Joint 

Differentiating  Eq.  B.9  once  with  respect  to  state  variables  gives 
the  following  nonzero  third  derivatives: 

3 

3  4> 

_ —  =  R2 

3 

30 

i 
3 

3  <f> 

_ El  =  -  R4 

3 

36 

j 


3 

3  (J) 

_ Z  =  -  R1 

3 

39 

i 
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3 

3  <)> 

_ y 

3 

30 

j 


=  R3 


B. 3.3.2  Translational  Joint 

Differentiating  Eq.  B.10  once  with  respect  to  state  variables  gives 
the  following  nonzero  third  derivatives: 


3 

3  <j> 


n 


3x  30 
i  i 


=  -  T1 


3 

3  <|> 


n 


3y  30 
i  i 


=  -  T2 


3 

3  <j) 


n 


3x  30 

j  i 


=  T1 


3 

3  <|> 


n 


3y  30 

j  i 


T2 


3 

3  <f) 


n 


30  30 
i  i 


=  -  T6 
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B.4  Cross  Derivatives  With  Respect 
to  Design  and  State  Variables 


B.4.1  Second-Order  Cross  Partial  Derivatives 

B.  4.1.1  Revolute  Joint 

Differentiating  Eqs.  B.7  once  with  respect  to  design  variables  the 
following  nonzero  second-order  cross  partial  derivatives  result: 


2 

3  <| > 


39  3b 
i  k 


2 

3  <j> 


39  3b 
j  k 


.  2 
3  <(> 


39  3b 
i  k 


2 

3  <}> 


39  3b 
j  k 


2 

3  $ 


9 


39  3b 
i  k 


2 

3  <j> 


9 


=  -  RB2 


=  -  RBI 


=  -  RB4 


=  -  RB3 


=  -  TB5 


=  TB5 


(B  .13) 


39  3b 
j  k 
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B.4»l«2  Translational  Joint 

Differentiating  Eq.  B.8  with  respect  to  design  the  following 
nonzero  second-order  cross  partial  derivatives  result: 
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B.4.2  Third-Order  Cross  Partial  Derivatives 


B.4.2. 1  Revolute  Joint 

Differentiating  the  second-order  derivatives  in  Eq.  B.9  once  with 
respect  to  design  gives  the  following  nonzero  third-order  cross 
partial  derivatives: 
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